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ABSTRACT 



A method is developed for digital simulation of linear time- 
invariant dynamic systems with lumped parameters and time delays. 
Ordinarily, such systems can be described by a linear matrix differential- 
difference equation, which can be transformed to an in finite -dimensional 
difference equation whose solution is obtained in a recursive way. 

As the present method depends on the accuracy of evaluation of the 
matrix exponential, a simple comnutationaj, procedure based on the 
truncation of the infinite series for e is described. 

In addition, an algorithm is given that ensures that the 
transient state of an unforced linear time -invariant dynamic system with 
zero time delay is calculated to a specified accuracy. 

Several sample problems are included. 
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CHAPTER 1 
INTRODUCTION 

1-1 Description of the problem 

This report presents a method for the simulation of linear time- 
invariant dynamic systems with lumped parameters and time delays. 

In many industrial processes one often encounters a type of time 
delay called "transportation lag". This kind of delay is generated when 
process materials move from one point in a process to another point 
without any appreciable change talcing place In the properties or 
characteristics of the process materials. Such delays may be caused by the 
flow of fluids through pipes, or by the motion of webs or filaments. 
Systems such as distillation columns and long heat exchangers are 
characterized by a multitude of small lags, which have an effect somewhat 
similar to that of time delays. The effects are not identical; however, 
some insight may be gained by using time delays models. The control of 
composition in a chemical reactor has been selected as a typical problem 
and this is depicted in section 5-2. 

Models having delays often arise in the study of systems with a 
mixture of lumped and distributed elements. An interesting form of 
topological representation suitable for such systems has been invented by 

Prof. H. M. Paynter at M.I.T. , and is called the bond graph. Rosenberg (17) 

it 
and Auslander (1) describe its use in modeling in some detail* 

Many other physical systems, such as electrical, mechanical and 



* Numbers in parenthesis refer to items in the bibliography. 



hydraulic transmission lines, and certain types of structural problems, 
are good examples of distributed systems which can be modeled using the 
delay operator. These systems often are analyzed as two-port chains, and 
usually the equations are slightly more involved than the type treated in 
this report. It is suggested that the reader interested in these kind of 
problems consult Koepcke (9) and Vaughn (20), as well as any standard text 
treating transmission phenomena. 

1-2 Formulation of the approach 

As an extension to the use of ordinary differential equations 
which arise when the future behavior of the system depends only upon its 
present state and not upon its past history, many systems that include 
time delays can be described by a linear matrix differential-difference 
equation. That is, the system is described by 

X(t) - I A. X(t - T ) + I D U(t - T ) , 
i-1 x X j-1 J J 

where X and U are the state and input vectors, respectively and T^ and T. 

are some fixed delay times. A. are a set of n x n matrices, and D are a 

set of n x r matrices. Techniques such as the direct method of lyapunov or 

laplace transforms can be used in the analysis of the equation. However, 

the use of these techniques frequently requires extensive computation, and 

for that reason they are not practical for hand analysis. At this step, 

designers and analysts are forced to rely on the digital computer as a 

computing aid. 

Because matrix manipulations are so convenient to implement on a 

digital computer, many existing dynamic systems programs are based on a 



matrix formulation of the problem, This convenience, together with the 
inherent elegance of the matrix approach, is helping to promote its 
acceptance among systems theorists. 

This report analyzes systems governed by the following differential- 
difference equation, for which it is desired to have a time sampled version 
of the state response: 

X(t) - A X(t) + B X(t - T) + D^Ct) + D^Ct - T) 
where 

T - time delay. 

X(t) - (n x 1) vector. It is called the state vector. 

U(t) » (r x 1) vector. It is the forcing signal or 

input vector, and it is assumed to be constant 
between samples. 
A, B ■ (n x n) constant coefficient matrices, 
D , D„ - (n x r) constant driving matrices. 

Koepcke (9) shows that the equivalent difference equation is (see 
section 3-1) 

00 

X(t + T) - I [*,(T) X(t - iNT) + A (T) U(t - INT)] 

i-0 

T 
where N - — , and 4> and A. are called plant transition matrices and 

control transition matrices respectively. 

The accuracy of evaluation of these sets of transition matrices 

depends upon the accuracy of evaluation of the matrix exponential. In 

section 2-3 a simple procedure based on the truncation of the infinite 

series of e (11,6), which guarantees a specified accuracy in the matrix 

exponential, is described. 



Also, a procedure is developed (21) to ensure that the calculated 
transient state of unforced linear time- invariant dynamic systems with 
aero time delay, is accurate to a specified tolerance. 

Several sample problems are presented to demonstrate the 
computation techniques. 

1-3 Application of results in dynamic simulation 

The two sets of simulators deduced throughout the development of 
this work, were tested on the time-shared IBM 7094 operated by Project MAC, 
and the entire operation, input and output, was carried out at an IBM 1050 
remote console typewriter, the algorithms will be part of the ENPORT 
Project which is being carried out at the mechanical engineering department 
under the direction of Professor Rosenberg. 

ENPORT is a digital computer program that accepts a bond graph 
description of a dynamic system and produces its time response. Work is 
being done on the theory of bond graphs, and a systematic graphical 
method has been developed for generating the state differential equations. 
ENPORT is organized in such a way that a broad class of nonlinear, active 
and passive, mixed energy- type systems can be handled. 

The wakelike nature of certain types of distributed systems make 
simulation by means of the digital computer, with its ability to exactly 
model the time delay operator, very natural. A simulation method based on 
delay-bond modeling has been developed by Aus lander (1) . 



CHAPTER 2 
DYNAMIC RESPONSE OF LINEAR TIME-INVARIANT SYSTEMS 

The analysis of many systems problems encountered In scientific 
and engineering investigations can be performed by either one of two 
major approaches. The essentially block diagram approach, Involves the 
determination of the transfer characteristics of the system components 
and the overall transfer characteristics. The second approach is based 
primarily upon the characterization of a system by a number of coupled 
first order differential equations which govern the behavior of the state 
variables. This technique is often implemented with the aid of a state 
variable diagram and is referred to as the state-variable approach. 

2-1 System Characterization by State Variables 

From the point of view of system analysis it is convenient to 
classify the variables which characterize or are associated with any 
system into (1) input, or forcing signals, Ui, which in essence represent 
the stimuli generated by systems other than the one under investigation 
and which influence the system behavior; (2) output, or response, 
variables Yi, which describe those aspects of system behavior that are 
of interest to the investigator; and (3) state variables Xi, which 
characterize the dynamic behavior of the system under investigation. 

One way of defining state variables is by making use of the state 
variable diagram. A state variable diagram is made up of Integrators, 
coefficients and summing devices. It describes the relationships among 
the state variables and provide physical interpretations of them. The 



outputs of the integrators denote the state variables. 

For continuous-time systems the state variable diagram is the 
same as the analog-computer simulation diagram. The state variable 
diagram may be derived from the overall transfer function of the system 
in three different ways (1) direct programming, (2) parallel programming, 
and (3) iterative programming. These methods are later ilustrated in the 
chapter corresponding to the solution to sample problems. Further 
information can be obtained from Tou (19), Schwarz and Friedland (18) and 
Ogata (15). 

2-2 Digital Solution of the Matrix Differential Equation 

A linear time-invariant system or process can be described by a set 
of first order linear differential equations with constant coefficients, 
which may be expressed in matrix form as 

X(t) - A X(t) + D U(t) (2.1) 

where 

A is the coefficient matrix 
D is the driving matrix 
X is the state variable vector 
U is the state forcing signal vector 
By analogy to the scalar case, the solution of eq. (2.1) is 

XCT) . e A < T " VxCt,,) + JV< T " *> B U( T ) dT (2.2) 

o 

with the initial conditions given by X(t ). 

For simplicity let t - 0, and let us define 



AT 

<KT) - e (2.3) 



as the transition matrix of the process. An equivalent name is the matrix 
exponential. 

Therefore eq. (2.2) can be reduced to 

X(T) - *(T) X(0) + *(T) e" AT D U(t) dT (2.4) 

Jo 

If X is small compared to the shortest period of interest in U(t), 
U(t) may be approximated over the region by U(0). 
Then eq. (2.4) becomes 

X(T) - *(T) X(0) + $(T)[ f e" AT dxj D U(0) (2.5) 



fT 

I e' AT dT - A" 1 [l - *(-T)J (2.6) 



-At 
By integration of the series of e 

f 

'o 

Thus 

X(T) m *(T) X(0) + *(T) A _1 [l - *(-T)] D U(0) (2.7) 
Let us define 

ACT) - te AT A" 1 - e AT A -1 e" AT ] D (2.9) 

as the control transition matrix. 

-AT 
From the series definition of e , It is observed that 

A -l ,-AT . e -AI A -l 

Therefore, eq. (2.9) becomes 

A(T) - [e AT A" 1 - e AT e" AT A" 1 ] D 

A(T) - [e AT A" 1 - A" 1 ] D 
or 

A(T) - [(e AX - I) A~ X ] D (2.10) 



Thus eq. (2.8) can finally be written as 

X(T) - e AT X(0) + [(e AT - I) A* 1 ] D U(Q) , (2.11) 



or 



X(T) - $(T) X(0) + A(T) U(0) (2.12) 

In general eq. (2.12) can be expressed as 

X(SflT) - *(T) X(KT) + A(T) U(KT) (2.13) 

which indicates that the state vector of the process after a particular 
interval depends upon the previous vector and also depends upon the 
forcing vector evaluated at the previous time. 

There are several methods available for computing the closed form 

AT 
expression for e , either as a special case of the study of the functions 

of a matrix or by a purely algebraic method based on the Laplace Transform. 

It is suggested, for those interested in these schemes, that they consult 

Ogata (15), Zadeh and Desoer (23), or Bellman (2). 

2-3 Digital Evaluation of the Matrix Exponential 

At 
e is given by 

e AT - e B - I + B + |< f,) + §( fj- ) + (2.14) 

note that each term in parenthesis is equal to the previous term. This 
provides a convenient recursion scheme. 

To ensure a reasonable truncation of the series, it is necessary 
to judge the convergence of the series. The norm of a matrix A is a real, 
non-negative number denoted by II A II , that gives a measure of the size of 
the matrix elements. 

Let 



*(t) - e AT « M + R 
inhere M is the truncated matrix which is an approximation of e (see 
reference 11) 

K /4 vi 



M " I ^f" (2 * 15) 

i-0 Xl 



and R is the remainder matrix 

ML' 1 

i-K+1 



R- J ^ (2.16) 



At 
If each element in the matrix e is required with an accuracy 

of "d" significant digits, then 

|r 1; ,|< llf^l (2.17) 

where r. . and m. , are elements of matrices R and M respectively. 
Let us define the norm of matrix A as: 

IaII - mindnax [J|a.,|] f max [[|a,,|]} (2.18) 

i J 3 J i 3 



(2.19) 
(2.20) 

(2.21) 
(2.22) 



and 



For 


this 


norm, we have 
II A Bll; 


cMIIbI 
|*|aI 








IIaIU-IIbIU i 


IaII + 


1b 


Then 


., it 


follows that 








U 


i«llv cat)*!! 1 


f HaII 1 


i 

T 



if the same norm is applied to the remainder matrix R. 
Upon expansion of eq. (2.22) 
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, i || A ||K+1 K+l || A |K+2 K+2 

I *,,,!< II A| j. II Al j , . 

1J = (K+l)! + (K+2)! + U,ZJ ' 

and, calling c the ratio of the second term to the first 



that 



UP 2 t k+2 



(K+2) 1 H aHt 
B . |K+1 K+l K+2 

B Al T 

(K+l)! 



(2.24) 



Therefore 

I UIIt . r (2.25) 

K - C 

Making the substitution of eq. (2.19) into eq.(2.23), it follows 



J K! K+l K! (K+2) (K+l) 

+ I1 <at) k 1I I (at^ || _ ( 26) 

+ K! (K+3)(K+2)(K+1) + * U ^° J 



Ir..|< || (At) K || [1 JAt I | At II II At 1 , 
3 K! ^ K+l K+2 K+l 

| At 1 1 At | |atJ 1 (2 2?) 

+ K+3 K+2 K+l r J \*-*n 



Thus 



|r t ,|< l(AT? K [ r[AT_L 1+ IAJ_J + 
1] K! [ K+l l L K+2 



||_At — HILAt — II , I (2.28) 

K+3 K+2 ' ' u " J 



Now, because any factor of the form ' Z TT for a>2 is always less 
than €, by eq. (2.24), then 
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r„l< «IAl£l(iy {1 + . + c \ t \ e \ _, 1 (2 . 29) 



KI { K+l 
If e<l, eq. (2.29) takes the form 

|r ± J<lkAl23f|Aj 1 

13 K! K+l * 1- e 



(2.30) 



This equation is suggested hy Everling (6) as the upper hound in 
the remainder matrix R. 

In order to Initialize the procedure, a certain K has to be 
chosen, but this K cannot be arbitrary, because it may happen that e>l, 
and relation (2.30) would not hold any more. 

This situation can be solved using eq. (2.25); thus 

K > IMI 

= e 
In order to ensure that e<L/2, the initial condition for K should 
be 

K> 2 II At II (2.31) 

However, it is possible that || Axil be less than 1/2; then K would 
be less than one. So, In order to avoid this possibility, an initial value 
of K can be obtained from 

K - max [ 2 1 At II , 2 ] (2.32) 

At this point, Everling (6) suggests that K be Incremented by 
half of its initial value, In the course of iteration. 

Although the matrix series approach for the evaluation of the 

transition matrix is suitable for digital computation, the disadvantage 

At 
stems from the convergence requirements for the series e , so it would 

be desirable to speed the computation. 

This can be done recalling the basic relationship 
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At j At /a 
e ■ e 



a 

(2.33) 



where a is chosen using the following expression 



a = l\ max la.. It (2.34) 

i,j U 3I J 

where 6 is the smallest integer allowed. 

The idea is to compute e T , because the norm of At /a is smaller 
than At, and the series will converge faster. Once the addition of the 
corresponding elements in the matrix terms of the infinite series is done, 
all that is required is to raise the result to the power a. The last step 
involves very few matrix multiplications, because a is a power of 2; for 
example, if a - 32 only 5 matrix multiplications are performed at the end 
of the computation. 

The steps presented in this section are summarized in a flow 
diagram in chapter 4. 

2-4 Error bounds in the transient response 

AT 

Although the matrix e can be obtained within prescribed 
accuracy, the truncation error of the matrix series, and the roundoff 
error do propagate in the state vector with increasing time. 

It is desirable, therefore to derive recursion relations which 
bound the propagated error due to these sources. Whitney (21) suggests 
one method. 

The homogeneous case of eq. (2.13) is 

X(K+T T) - *(T) X(K T) (2.35) 

If eq. (2.15) is used in place of *(T), the numerical calculation 
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reads 

X^(K+T T) « M X^(K T) (2.36) 

vhere X^(K+1 T) is the perturbed state vector obtained from numerical 
calculation. 

The propagated error at time (K+l) T due to the approximate M 
is 

E(K+T T) - X(K+T T) - X^(K+T T) (2.37) 

Rewriting eq . (2.35) and substracting eq. (2.36) from it yields 

X(K+T T) - X^(K+T T) - [MfR][X*(K T) + E(K T) ] 

- M X^(K T) , (2.38) 
or 

E(K+T T) - [M + R]E(K T) + R X^(K T) (2.39) 

From eq. (2.17) 

Ir^HTVyl (2.17) 

We can define 

** = |r ±j | I (2.40) 

where I is a matrix each of whose elements is unity. Replacing R with R^ 
in (2.39), we obtain the running error bound for E(K+1 T) , that is 

E(K+T T) - [M + R] E(K T) + R* X* (K T) (2.41) 

The computation may be initialized assuming E_(0) is zero. 
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CHAPTER 3 

DYNAMIC RESPONSE OF LINEAR TIME- INVARIANT SYSTEMS 
WITH LUMPED PARAMETERS AND TIME DELAYS 

It has been found that many industrial processes in which 
transportation lags are common can be described by a system of 
differential-difference equations. The chemical process industry offers 
many examples. 

This chapter analyzes the special case of a system subject to 
one delay, and a technique suitable for digital computation is described. 
The derivation follows a criterion developed by Koepcke (9). 

3-1 Digital solution of the matrix differential-difference equation 

Consider a dynamic system which is governed by the following 
differential-difference equation 

X(t) - A X(t) + B X(t - T) + D^Ct) + D^Ct - T) (3.1) 

X(t) - (n x 1) vector, referred to as the state vector; 
U(t) - (r x 1) input vector, assumed constant between samples; 
i.e. U(t) - U(t k ) for t^t^tj^j 

A, B - (n x n) constant coefficient matrices; and 
D., D„ - (n x r) constant driving matrices 

Let us consider first the homogeneous part of eq. (3*1); that is 

X(t) - A X(t) + B X(t - T) (3.2) 

Taking the laplace transform of eq. (3.2), 



where 
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SX(S) - X(0) - (A + B e ai ) X(S) (3.3) 



X(S) - [SI - (A + B e" 81 ]" 1 X(0) (3.4) 

Defining Z=e , than 

X(S) - [SI - (A + B Z)]" 1 X(0) , (3.5) 



X(S) - | [I - (A + BZ)/S] -1 X(Q) . (3,6) 

Let W - [I - R]" 1 , where R - A ^ BZ j then 

W-I + R+R 2 + R 3 + R 4 + . (3.7) 

Therefore, one should choose an "S" large enough to ensure that 
eq. (3.7) la valid. 
Thus 

X(S) ,A [1+ A^Z + (A + BZ} 2 + <A + BZ? 3 + 

s s s 2 S 3 

+ < A + BZ > + J X(0) (3.8) 

S* 

Recall the facts that 

(A + BZ) 2 - A 2 + A(BZ) + (BZ)A + (BZ) 2 

(A + BZ) 3 « A 3 + A 2 (BZ) + A(BZ)A + A(BzJ + (BZ)A 2 + 

+ (BZ)A(BZ) + (BZ) 2 A + (BZ) 3 

(A + BZ) 4 - A 4 + A 3 (BZ) + A 2 (BZ)A + A 2 (BZ) 2 + A(BZ)A 2 + 

+ A(BZ)A(BZ) + A(BZ) 2 A + A(BZ) 3 + (BZ)A 3 + 

+ (BZ)A 2 (BZ) + (BZ)A(BZ)A + (BZ)A(BZ) 2 + 

+ (BZ) 2 A 2 + (BZ) 2 A(BZ) + (BZ) 3 A + (BZ) 4 



etc. 



Then, arranging by terns of equal delay, 
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I A A 2 A 3 A 4 

X(S) - [± + £+£- + £- + £- + 1 X(0) + 

b s z S J S* S^ 

. r BZ . A(BZ) + (BZ)A A 2 (BZ) + A(BZ)A + (BZ)A 2 
1 2 3 4 

+ A 3 CBZ) + A 2 (B Z)A + A(BZ) A 2 + (BZ)A 3 + _ } ^ (Q) + 
S 5 

, r (BZ) A(BZ) 2 +(BZ)A(BZ) + (BZ) 2 A A 2 (BZ) 2 + A(BZ)A(BZ) 
[ S 3 S 4 s 5 

+ A(BZ) 2 A + (BZ)A 2 (BZ) + (BZ)A(BZ)A + (BZ)V | 
S 5 

+ -- 1 x(o) + 

. r (BZ) 3 , A(BZ) 3 4- (BZ)A(BZ) 2 + (BZ) 2 A(BZ) + (BZ) 3 A 
1 4 5 

+ — ] X(0) + 

+ [^1S2- + ] x(0) + 



Now, because 



We have 



+ (3.9) 



ZX(t) - X(t - T) (Z = e" ST ) , (3.10) 



ZX(0) - X(-T) , 

Z 2 X(0) « X(-2T) ( 

Z 3 X(0) - X(-3T) , and so forth. 
Therefore, X(S) can be arranged In the following way. 

X(S) - *o(S)X(0) + *i(S)X(-T) + $ 2 (S)X(-2T) + * 3 (S)X(-3T) + 

+ **(S)X(-4T) + (3.11) 



where 
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2 3 4 

$oCS) - - + ^T + ^T + \ + ^T+ (3.12) 

S S" 1 S J S 4 S 3 

* ^n B j_ AB + BA _,_ A 2 B + ABA + BA 2 ^ 
*i (S) - -r + r— + 7 + 

j A 3 B + A 2 BA + ABA 2 + BA 3 ; (3.13) 

s 5 

* /o\ B 2 ^ AB 2 + BAB + B 2 A ^ A 2 B 2 +ABAB + AB 2 A + BA 2 B , 
<& 2 (S) - -r + r + c + 

S J S* S 

+ BABA + B 2 A 2 + (3t u) 

s 5 

«,<S) - 4 + Afl3 + ^ t ^^ + B ' A + (3.15) 

iL 

s 5 

Rearranging terms, It follows that 

T A a 2 a 3 a* 

MS) - i + A r + A r + A r + A r + (3.17) 

S S^ S S S 

. ,_ s B A AB + BA A A(AB + BA) + BA 2 . 
*! (S) - — + — + -* f + 

+ A[A(;AB + BA) + BA 2 ] + BA 3 + (3#18) 

s 5 

$2 ( S ) . *L + AB 2 + B(AB + BA> + A[AB 2 + B(AB + BA) ] + 

s 3 s 4 S 5 

+ B[A(AB + BA? + BA 2 ] + (3>19) 

s 5 

* ,(3) - 4 + AB3 + B t AB2 s +MAB + BA)] + » (3.20) 

B 4 
*„(S) - ^r + (3.21) 

Let us try to find a relationship among the coefficients. With 
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this idea in mind, we shall form the array shown in figure 3.1: 
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It is seen that the correlation among the elements (call any 



element by C . . ) Is 



s" 1 S~ 2 S"" 3 S -4 S -5 



r" 



ol 



-11 



: 12" 



u o2 



-** 



"22 



_* 



"ol 



13 



23 



-^ 



"33 



o4 



"14 



"24 



"^ 



34 
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where the arrows Indicate the inmediate dependance; i,e. t C^ 2 depends on 



C _ and C 11( etc. 
ol 11 



From a careful study of the array In fig. 3.1, it Is found that 

(3.22) 



C. . AC... B C, - . , 



S* s* +1 



J+l 



where "1" is the subindex denoting row and "j" is the subindex denoting 
column. 

The following conditions should be added, in order to initialize 
a computational procedure 



°-i.j ' ° & 



C I 

o.o ■ 



C 4 - i>0 
i,o 



(3.23) 
(3.24) 
(3.25) 



The inverse laplace transform of eq. (3.22) yields (note: 
L[t n /n!] - 1/S n+1 ) 
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[c i,j ] oir - [ AC i. 3 -ii ir + ' Bc i-i (j -i] "TT (3 - 26) 



Therefore 
C 



i A " tAC i i-l ] TJ (j " 1)! + t BC i-l ^JjilblU -C3- 27 ) 
,J ,J i jl T J-1 ,J jl T J-1 



[AC Jt + [BC ]t 

c - iU-i - 1 A ll ■ + ■ . (3.28) 

^pJ 3 

Changing j for j+1, eq. (3.28) Cakes the final form 

[AT]C + [BT]C, , 
C i,j+1 J+1 U "" ; 

Actually eq. (3.29) gives all coefficients without any need to 

multiply them by r^ 

jl * 

This is because T has been associated with matrix A and B f and in 

order to compute any C. , -, the initial conditions given by eqs. (3.23), 
l,j+l 

(3.24) and (3.25) have to be considered. 

The computation of the C, . . is done in a recursive way, as 
I, j+1 

given by eq, (3*29). Once they are computed, they may be substituted in 
the inverse laplace transformation of eqs. (3.17), (3.18), etc., so that 
*o(t), *i(t), *2(t), ... can be generated. The last set of matrices are 
called "plant transition matrices". 

Returning to eq. (3.11), if e C is multiplied into both aides, 



then 



e CS X(S) - $ (S)e tS X(0) + ^ (S)e tS X(-T) + $ 2 (S)e tS X(-2T) + 

+ $ 3 (S)e tS X(-3T) + , (3.30) 
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e tS X(S) - * (S)X(t) + *i(S)X(t - T) + * 2 (S)X(t - 2T) + 

+ $ 3 (S)X(t - 3T) + . (3.31) 

Taking the inverse lapXace transform of eq. (3.31), it turns out 
to be 

X(t + T) - *o(T)X(t) + *i(T)X(t - T) + * 2 (T)X(t - 2T) + 

+ MT)X(t - 3T) + , (3.32) 

or 

00 

X(t + X) - I * t (T)X(t - IT) (3.33) 

1-0 
This is the sampled version of the homogeneous part of the 
differential-difference equation. 

Now, let us consider the addition of an input vector or forcing 
signal. 

In chapter 2, section 2-2, it was found that the digital version 
of the time-invariant matrix differential equation adopted the form 

X(Rfl T) - *(T) X(KT) + A(T) U(KT) (3.34) 



where 



and 



*(T) - e 



, AT (3.35) 



A(T) _ (e AT - I) A" 1 D (3.36) 

Although it was not demonstrated, it can be shown that 



Am - y < AT ? J 
flW i (j+i)t 

j-0 


T D , 


MX)- I U$- 

j-0 


— 1 — T D 
j+1 



(3.37) 



(3.38) 
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If the terms -rrr and TD were absent, the series would be the well 
known matrix exponential, whose terms can be computed in a recursive way by 



C o>j - jfil (3.39) 

Therefore, eq. (3.38) is 

ACT) -j %i ik ° (3 - 40) 

By following the same line of reasoning, the control transition 
matrices in the case of the complete differential-difference equation can 
be written as 

CO oo 

4 i <T) " j^l.J J& D i + ^"i-LJ 3^ D 2 »- 41 » 

and the complete difference equation is 

00 

X(t + x) - I [*.(T) X(t - IT) + A.(T) U(t - IT)] (3.42) 

i-o 

In resume, the digital version of 

X(t) - A X(t) + B X(t - T) + D,U(t) + DJJ(t - T) 



is 



where 



X(t + T) - I [4.(T) X(t - iNT) + A (T) U(t - INT)] 

i-0 X 



T 



1 j-o i,J 



[AT] C, . + [BT] c, , , 

i.J+1 j+1 



C - I 

o,o 
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CHAPTER 4 
ALGORITHMS FOR DIGITAL COMPUTATION 

This chapter presents flowcharts for the algorithms of chapters 
2 and 3, from which the computer programs were derived. They accept as 
Input the coefficient matrices, the driving matrices, the initial state 
vector, and deterministic forcing vectors. As output, the computer will 
produce the state vector at the current sampling time and the set of 
transition matrices, if desired. 

Because these routines will eventually become part of Project 
ENPORT, they were designed to be used on the time-sharing system. However, 
they may be operated in the BATCH procedure without any difficulty, by 
modifying the input/output statements. 

The programs were written in the MAD language, and are listed in 
Appendix A. 
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4-1-1 TRANS 



Purpose: to compute the time response of linear time- invariant 
systems. 



Inputs: order of system (M - ) ; sampling time (T - ) ; final 
time (TF - ) ; number of input signals (R - ) ; the 
augmented A matrix and the initial state (x(l) - ). 

Outputs: the transition matrix; the current time; and the state 
of the system. 

Remarks: main program. Subroutines called by TRANS: EXPMAT, and 
DISTUR. 
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Read from the console 
order of system (M - ) 
sampling time (T - ) 
final time (TF - ) 



Is there any disturbing 
signal ? 



yeu 



Read from the console 
number of input signals 
(R - ) 



M - M + R 



Read from the console 
the augmented A matrix 
A(l,l) - ~, A(2,l) - - 
A(M,1) - — 



Read from the console 
initial state X(l) - - 



Save initial state 
XI(1) - X(l) 
XI(M-R) - x(M-R) 



TA - T 



~W 



R - 



TA is running time 
TRANS. Page 1 of 3 pages. 



i 



Is R + ? 



yes, 



Execute distur. (TA) 



XKM-R+1) - X(M-R+1) 



XI (M) - X(M) 



Initial error 

E(l) - 0. 

E(M) - 0. 
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E(I) is error vector 



TZ - T 



Get transition matrix 
Execute expraat. (T) 



Output to console 
Transition matrix 
EM(I.J) 



Save transition matrix 
EMP(I,J) - EM(I,J) 



Initialization 

PE(I) - 

Y(I) - 



TZ is time increment 



<*— © 



PE is new E 



Y is new X 
TRANS. Page 2 of 3 pages 
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Compute new state 
Y(I) - Y(I) + EMP(I,J)*X(J) 

I 



PE(I) - PE(I) + (EMP(I.J) + RIJ)*E(J) 
+ RIJ*X(J) 



RIJ is upper bound in 
remainder terms 



NORM - I |PE(I) | 
I 



Is NORM > 10™ 7 ? 



yes 



T - TA 

~T" 



Execute expmat. (T) 



PE(I) - 0. 
Y(I) - 0. 



Y(I) - Y(I) + EM(I,J)*XI(J) 
PE(I) - PE(I) + RIJ*XI(J) 



X(I) - Y(I) 
E(I) - PE(I) 



Output to the console 
TA, X(l) ... X(M-R) 



, Is TA < TF ? 

_3 



na 



TA - TA+TF 



Ayes 



Get new disturb, vector 
Execute distur. (TA) 



End of program 



J 



TRANS. Page 3 of 3 pages. 
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Entry to expmat. 



B(I,J) - A(I,J) 



B(I,J) - B(I,J)*T 



a * 2 > nu*|B(I,J)j 

3 is the smallest positive 

Integer 



B(I,J) - B(I,J)/a 



NORM- |B| - «in{max[i:|B(I,J)|], maxIElBCl.J)!]} 
I J J I 



Initial values of K 
K - uax[2*N0RM, 2] 



K is the number of 
terms of the series 



Increment 
IN - K/2 



K B 1 



1-0 i! 



T 



-© 



e - 



IBI 



K+2 



I 



EXPMAT. Page 1 of 2 pages. 
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© 
i 



TJT T 




b k r 

Kl 


IIbII i ] 






k+l l-e J 






v 








la RIJ < 10~ 7 |EM(I,J)| ? 


no 


K - K + IN 










yas 


u 






1 




U 


. - l 


© 














i r 






TERM - EM*EM 










u 








LL - LL + 1 










V 












Is LL>0 ? 


no 


EM - TERM 












y.s 


' 












Function return 




1 


r 






End of function 





EXPMAT. Pag. 2 of 2 pagaa. 






the forcing signal ve:tor at the current 



utine '"-alJ.ea ; v r ^A.-i.j, 



L',XP^VI 



cne matrix expoi'^r. I i.a i. , 



tkan:-;. 
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_ ' , i 



jf_ 

M - i< 



_j „„ f 



fK) 



1 

;;;t ion return | 



.» n 

:>f t unction : 



DISTUH., Page. ,, or A page. 
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4*2*1 TIMDEL 

Purpose: to compute the time response of linear systems with 
lumped parameters and time delays. 

Inputs: order of system (M - ) j sampling time (T - ) ; time 
delay (TD - ) ; final time (TF - ) ; number of input 
signals (R - ); the A matrix; the B matrix; the D^^ 
matrix; the D matrix; the initial state (X(l,l) - ). 

Outputs: the plant transition matrices, the control transition 
matrices if desired; the current time; and the state 
of the system. 

Remarks: main program. Subroutines called by TIMDEL: DELFOR, 
and PERTUR, 
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Read from the console 
order of system (M - ) 
sampling time (T - ) 
time delay (TD - ) 
final time (TF - ) 



Is there any disturbing 
signal? 



yes 



Read from the console 
number of Input signals 
(R - ) 



REL - TD/T 



Read from the console 
the A matrix. A(l,l) 
A<2,1) - — 



Read from the console 
the B matrix. B(l,l) - 
B(2,l) 



no 



r - o 



REL > l t Integer 



i" 



Is R ¥ ? 



y«« v 



Read from the console 

the D. matrix. D^l.l) - -- 

D 1 (2,l) - - 



Z 



TIMDEL. Page 1 of 3 pages. 
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©- 



©- 



1 



Read from the console 
the D 2 matrix. D (1,1) 
D 2 (2.1) - - 



Compute set of transition 

matrices 

Execute delfor. (T) 



Do you wish to have the 
transition matrices ? 



yes 



Output to the console 
plant transition matrices 
BM(L) , control transition 
matrices DELF(L) 



Bead from the console 
initial state X(l,l) - 



W is the number of 
transition matrices 
computed 



X is (M x W) 



TA - 



W-l 

X(TA+T) - I EM(I+1) X(TA - I*REL*T) 
1-0 ~ 

I 



- [is R + 7| 

_Z!5; 



Execute pertur. (TA) 



~^l 

X(TA+T) - X(TA+T) + J DELF(I+1) U(TA - I*REL*T) 
1-0 

3 



TA - TA + T I 



Output to the console 
TA, X(TAfT) 



KD 



TIMDEL. Page 2 of 3 pages. 



it\*~t- • u' j- 



:ogi 8.13 






:,)jU x: . s e p^ant trans, r 
il transition matrices, 



i:itn.ce-s ana t:\e 



<:. i:\ti mailed bv IIMDLI. 
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Entry to delfor 



A(I,J) - A(I t J)*T 



B(I,J) - B(I,J)*T 



I 



Is R 4 ?~| 



yea 



D 1 (I,J) - D 1 (I,J)*T 



I 



D 2 (I f J) - D 2 (I,J)*T 
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NORM- ||A|| - Mln max[z|A(I,J)|], n,ax[j:|A(I,J)|] 
I J J I 



Initial value of K 
K - max[2*N0RM t 2] 



Increment 
IN - K/2 



i - o 

— T 



C(-1,J) - t J>0 

I 



C(0,0) - I I 



G<2,-1) - 



C(I,0) - 0, l>0 



K ia the number of terms 
of the seriea e AT 



■O 



DELFOR. Page 1 of 2 pages. 
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£ 



CU.Jfl) - [A]C(I,J) + [B]C(I-1,J) 



J+l 



± 



EM(I+1) - I C(I,J) 
J«0 



T 



0(i,O - I ^^2i D L 



J- 1 



J+l 



L - LrKL 



I 



Is L>2 ? 



yes 



DELF(I+1) - G(1,I) + G(2,I-1) 



K - K+IN 



Is I i* ? 



yes 



NORM 
K+2 



I 



RIJ - 



A K |ir N0RM 1 1 
KI ||Lk+1 1-e J 

2 



Is RIJ < 10" 7 |EM(I)| ?] 



yes 



Is ||em(i+dj| < io" 7 ? 



yes 



RIJ is upper bound In 

AT 
remainder terms of e 



I * 1+1 



Function return 



End of function 



DELFOR. Page 2 of 2 pages. 
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PEKTLR 

. ompute tut forcing si^nai sector at trie current 
.: .. irie pracram tias to keep track of the past. 



jtine called bv TIMDl'I.. 



;::**£« . 
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CHAPTER 5 
SOLUTION TO SAMPLE PROBLEMS 

This chapter describes a set of sample problems which were 
selected because they represent typical applications of the two simulators. 
They are intended to show the use of the state variable diagram, and also 
to show the accuracy of the methods, 

5-1 Test problem for the simulation of dynamic systems without delay 



Example 5-1. Although this example may represent a great number 
of physical processes, it was selected purely from the mathematical point 
of view. The same problem was run by Liou (11). 
Given 

1 0~ 
X(t) - 1 X(t) (5.1) 

-.75 -2.75 -3 



and 



X(0) - 



2 
-2.5 

3.75 



(5.2) 



Obtain X(nT) using T - 0.1 Min. 

The reported solution by Liou and the one obtained by the 
simulator are 
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loadgo trans expmat distur 

W 1010. U 

EXECUTION. 

GIVE ORDER OF SYSTEM (M = ) 

SAMPLING TIME (T = ), FINAL TIME (TF = ) 

m-3, t«.l, tf=2.* 

IS THERE ANY DISTURBING SIGNAL 
no 

GIVE THE A MATRIX ( AC 1, 1 )=--, A( 2, 1 ) =-- ) 

aCl,l)=0.,l.,0.* 

a(2,l)=0.,0.,l.* 

3(3,1)=-. 75,-2. 75,-3.* 

GIVE INITIAL STATE CX< 1)= — ) 
x(l)=2.,-2.5,3.75* 

TERMS OF THE MATRIX EXPONENTIAL 



v(t) 


- A V(t) 




v(t) - 


X(l) 
X(2) 
X(3) 




A- 


1 ' 
1 
-.75 -2.75 -3 


V 


CO) - 


2 
-2.5 
3.75 





EM{ 


1, 1) = 




99988UE 00 








EMC 


1, 2) = 




995717E-01 








EM( 


1, 3) = 




U52513E-02 








EM( 


2, 1) = 


-. 


339385E-02 








EM( 


2, 2) = 




9874U0E 00 








EM( 


2, 3) = 




859963E-01 








EMC 


3, 1) = 


-. 


6^972E-01 








EM( 


3, 2) = 




23988<*E 00 








EM( 


3, 3) « 




729i*51E 00 








TIME 


X(l) 




X(2) 




X(3) 




.10 


.176781E 


01 


-.215290E 


01 


.320616E 


01 


.20 


.15677UE 


01 


-.18561UE 


01 


.27U16E 


01 


.30 


.139515E 


01 


-.160242E 


01 


.234368E 


01 


.«*0 


.12U603E 


01 


-.1385U8E 


01 


.200401E 


01 


.50 


.111700E 


01 


-.119997E 


01 


.171382E 


01 


.60 


.100515E 


01 


-.104131E 


01 


.1U6596E 


01 


.70 


•907978E 


00 


-.905571E 


00 


. 125U31E 


01 


.80 


.823379E 


00 


-.789U13E 


00 


.107362E 


01 


.90 


.7U9538E 


00 


-.68996UE 


00 


.919U18E 


00 


1.00 


.68I+911E 


00 


-. 604775E 


00 


.787838E 


00 


1.10 


.628178E 


00 


-.531753E 


00 


.675590E 


00 


1.20 


.578215E 


00 


-. U69107E 


00 


.579853E 


00 


1.30 


.53U062E 


00 


-.M5312E 


00 


.1+98212E 


00 


l.UQ 


.i*9U901E 


00 


-.36906UE 


00 


.U28602E 


00 


1.50 


.46003«iE 


00 


-.329250E 


00 


.369257E 


00 


1.60 


. 423868E 


00 


-.29i*921E 


00 


.318666E 


00 


1.70 


.400894E 


00 


-.265268E 


00 


.275537E 


00 


1.80 


.375681E 


00 


-.239602E 


00 


.238768E 


00 


1.90 


.352861E 


00 


-.217331+E 


00 


.207H5E 


00 


2.00 


.332118E 


00 


-.197965E 


00 


.180676E 


00 


2.10 


.313185E 


00 


-.181068E 


00 


.157862E 


00 


END OF EXECUTION 












TO CONTINUE, 


GO TO THE TOP OF A NEW PAGE 








AND PRINT AN 


ASTERISK 













Figure 5.1 Console transaction for example 5.1 
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A liquid stream enters tank 1 (figure 5,3) at a volumetric flow 

3 
rate F cfm and contains reactant A at a concentration of C q moles A/ ft . 

Reactant A decomposes In the tanks according to the irreversible chemical 
reaction. 

A 1- B 

The reaction is first order and proceeds at a rate 
r - k c 

where 

3 
r ■ moles A decomposing/ (ft )(time) 

3 
c - concentration of A t moles A/ft 

k - velocity constant, a function of temperature 

The reaction is to be carried out in a series of two stirred 

tanks. The tanks are maintained at different temperatures. The temperature 

in tank 2 is to be greater than the temperature in tank l t with the result 

that k , the velocity constant in tank 2, is greater than in tank 1, V.^, 

Changes in physical properties due to chemical reaction are neglected. 

The purpose of the control system is to maintain c„, the 

concentration of A leaving tank 2, at some desired value In spite of 

variation In inlet concentration c . This will be accomplished by adding a 

o 

stream of pure A to tank 1 through a control valve. 

Further assumptions are that the control valve and the measuring 
element have no dynamics, and that the controller exert proportional action 
on the process, 

A portion of the liquid leaving tank 2 is continuously withdrawn 
through a sample line. The measuring element is remotely located from the 
process, because rigid ambient conditions must be maintained for accurate 
concentration measurements. The sample line can be represented by a 
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transportation lag* 



pure A 



i 




control la r 



Vr c i 



V » T o. c o 



-±r 



ff 



product 
■t re an 



composition 
measuring 



sample 
stream 



heating coils 
Figure 5.3 
Control of a stirred- tank chemical reactor 

The following data is assumed to spply to the system 
Molecular weight of A - 100 lb/lb mo la 
P A - 0.8 lb mole/ft 3 

C - 0.1 lb mole A/ ft 3 

OS 

F - 100 cfm 
m - 1.0 lb mole/min 



l^ - 1/6 min 



-1 



k 2 - 2/3 min' 



V - 300 ft J 



-1 



Valve sensitivity k - 1/6 cfm/psi 

Measuring device sensitivity 

k - 100 in. pen travel/ (lb mole/ft 3 ) 
m 

Time delay in sample line - T 

The overall block diagram which the authors propose is 
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•Q- 



S + 1 



2S + 1 



-TS 



Figure 5.4 

Block diagram for a chemical reactor 

control system 

It Is assumed that the inlet concentration c does not change 

o 

with time. 

As was discussed in chspter 2, the state variable diagram can 
be obtained in three ways. Direct programming will be used in this case. 
With this purpose, let it be called c A the input to the lag term ■— 
and c its output in figure 5.4, then 



c L_ 

c B 2S+1 • 



(5.3) 



c g ' 5 1+.5 S 



Eq. (5.4) can be written as 

c - .5 S" 1 



(5.4) 



(5.5) 



where 



1+.5S 



-1 



Transposing 



-1 



■b " C B - ' 5 S «b 



(5.6) 



(5.7) 
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By a similar procadurs 



1-JL 

c A ftfl 



(5.8) 



C B S 



-1 



C A 1+S" 1 



<B " S " 1 \ 



vhsra 



-1 



Tranapoalng 



1+8 



E - c. - S" 1 c. 
• A A 



(5.9) 



(5.10) 



(5.11) 



(5.12) 



Tha 



■riabla diagran follows from aqs. (5.5), (5.7) and 



aqa. (5.10), (5.12), and la shown In flgura 5.5. 

^(0) 

X 2 




Plgura 5.5 

Stata varlabla diagram for a chawlcsl raactor 

control ays con 

Iha notation in flgurs 5.5 haa boon changsd ■ lightly. Thla la In 
ordar to follow tha aano synbolisn glvan In tha pravioua chaptara. 

In flgura 5*5 tha atata varlablaa aro X, and Xj. Tha dlffarantlal- 
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difference equations for the state variables are readly obtained by 
Inspection of the diagram. That is, 

X x - -.5 X x + .5 X 2 



X 2 «-X 2 +KU-K X x (t - T) 
Therefore the matrix differential-difference equation is 



(5.13) 
(5.14) 



X(t) - 



-.5 


.5^ 


x(t) + 








x(t - 


- T) + 





Q 


-1 




-K 









K 



U(t) (5.15) 



where 



X(t) - 



x x (t) 
x 2 (t) 



from this equation, it is seen that the coefficient matrices and 
driving matrices are 

-.5 .5 

-1 



A - 



(5.16) 













B - 










_ " K 








D l" 




"o" 




D 2- 










(5.17) 



(5.18) 



(5.19) 

Five numerical examples were run using this system. These are 
summarized as follows. 

Example 5.2.1. Overall forward gain K- 5,24. We assume a time delay 
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equal to zero* A unit step is the input and all initial conditions are 
zero. The matrix differential equation is 



where 



V(t) - 



-.5 .5 

-5.24 -1 5.24 





V(t) 



(5.20) 



V(t) 



x x (t) 
x 2 (t) 

U(t) 



(5.21) 



Example 5.2.2. Overall forward gain K - 5.24, and time delay - 
.5 Min.. Same conditions of the state were taken. The matrix differential- 
difference equation is 



X(t) - 



-.5 .5 
-1 



X(t) + 




-5.24 






X(t - .5) + 



U(t - .5) 




5.24 



U(t) + 



(5.22) 



Example 5.2.3. Overall forward gain K - 1.85. Time delay is zero. 
The remaining conditions are the same. The state equation is 



V(t) 



-.5 .5 

-1.85 -1 1.85 





V(t) 



(5.23) 



Example 5,2.4, Overall forward gain K - 1.85. Time delay » ,5 min. 
Unit step and zero initial conditions are assumed. The state equation is 
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X(t) - 



-.5 .5 
-1 



X(t) + 




-1.85 





X(t - .5) + 



U(t - .5) 




1.85 



U(t) + 



(5.24) 



Example 5.2.5. This ia the same as example 5.2.4, with the 
exception of the time delay, which is taken equal to 1 rain*. The state 
equation is 



£U) - 



-.5 


.5" 


x(t) + 





o" 


£<* ■ 


- 1) + 


" 





-1 




-1.85 









1.85 



U(t) 



U(t - 1) 



(5.25) 



All five examples with the input/output information and the 
response curves, are shown in figures 5.6 to 5.15. 

The interested reader should compare the responses of the three 
cases with delay with those given by Coughanowr and Koppel on page 467 of 
reference (4) . 

5-3 Test problem for the simulation of dynamic systems with delays 

The eighth example was run in order to check the accuracy of 
evaluation of the set of transition matrices. This example is discussed by 
Koepcke (9). 

The problem is described as an unstable process which is governed 
by 



54 



GIVE ORDER OF SYSTEM (M = ) 

SAMPLING TIME (T = ), FINAL TIME (TF = )| 

m«2,t-.5,.tf»ll.* | 

IS THERE ANY DISTURBING SIGNAL j 

yes | 

GIVE NUMBER OF INPUT SIGNALS (R = ) 

r = l* 

GIVE THE A MATRIX ( A( 1, 1 ) =--, A( 2, 1 ) =-- ) 
3(1,1)=-. 5,. 5,0.* 
a<2,l)=-5.2t* / -l.,5.2i+* 

a<3,l)=Q.,0.,0.* 

GIVE INITIAL STATE (X(l)=--) 
x(l)=0.,0.* 

TERMS OF THE MATRIX EXPONENTIAL 



V(t) - A V(t) 

X(l) 
V(t) » X(2) 
U(l) 



EMC 


1, 1) = 




556076E 00 




EMC 


1, 2) = 




154089E 00 




EMC 


1, 3) = 




243387E 00 




EMC 


2, 1) = 


- 


161485E 01 




EMC 


2, 2) = 




401987E 00 




EM C 


2, 3) = 




185824E 01 




EM( 


3, 1) = 




O000O0E 00 




EMC 


3, 2) = 




000000E 00 




EMC 


3, 3) = 


1 


000000E 00 




TIME 


Xtl) 




XC2) 




.50 


.243387E 


00 


.185824E 


01 


1.00 


.665063E 


00 


.221219E 


01 


1.50 


. 9 5 4 8 7 E 


00 


.167353E 


01 


2.00 


.1031S1E 


01 


.990269E 


00 


2.50 


.969739E 


00 


.590102E 


00 


3.00 


.873563E 


00 


. 5 2946 8 F 


00 


3.50 


.810740E 


00 


.660403E 


00 


4.00 


. 795981E 


00 


.814488E 


00 


(+.50 


.811516E 


00 


.900262E 


00 


5.00 


•833372E 


00 


.909654E 


00 


5.50 


.846973E 


00 


. 8 78 136 E 


00 


6.00 


•S49679E 


00 


.843503E 


00 


6.50 


.845848E 


00 


.825210E 


00 


7.00 


. 8I+0898E 


00 


.824044E 


00 


7.50 


.837966E 


00 


.831568E 


00 


8.00 


.837495E 


00 


.839327E 


00 


8.50 


.838429E 


00 


.843206E 


00 


9. 00 


.839546E 


00 


.843258E 


00 


9.50 


.840175E 


00 


.841475E 


00 


10.00 


.840250E 


00 


.839743E 


00 


10. 50 


. 84002 5E 


00 


.838325E 


00 


11.00 


. 839 7 74 E 


00 


. S33960E 


00 


END OF EXECUTION 








TO CONTINUE, 


GO TO THE TOP OF A NEW PAGE 




AND PRI NT AN 


ASTERISK 









-.5 .5 

-5.24 -1 5.24 





XCO) - 10 



Figure 5.6 Console transaction for example 5.2.1 
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1 
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loadgo timdel delfor pertur 

W 10 3 9. U 

EXECUTION. 

GIVE ORDER OF SYSTEM (M = ) 

DESIRED SAMPLING TIME (T = ) 

TIME DELAY ( TD = ), FINAL TIME CTF = ) 

m=2, t = .5, td=.5, tf-15.* 

IS THERE ANY DISTURBING SIGNAL 
yes 

NUMBER OF INPUT SIGNALS (R = ) 
r = l* 

GIVE THE A MATRIX ( A( 1, 1)=--, A( 2, 1 )=-- ) 

8(1,1)=-. 5,. 5* 

a(2,l)=Q.,-l.* 

GIVE THE B MATRIX ( B ( 1, 1 ) =--, 3 ( 2 , 1 ) =-- ) 
bCl,l)=0., 0.* 
b(2,l) — 5.2U,0.* 

GIVE THE Dl MATRIX ( Dl C 1, 1 ) =--, Dl ( 2, 1 ) =" ) 

dlC 1, 1)=0.* 

dlC2,l)=5.2i+* 



GIVE THE D2 MATRIX (D2C1 

d2(l,l)=0.* 

d2(2,l)=0.* 

DO YOU WISH TO HAVE THE 

yes 

TRANSFER MATRIX 

EMC 1, 1) = 

EMC 1, 2) = 

EMC 2, 1) = 

EMC 2, 2) = 

TRANSFER MATRI X 

EMC 1, 1) = 

EMC 1, 2) = 

EMC 1, 1) = 

EMC 2, 2) = 

TRANSFER MATRI X 

EMC 1, 1) = 

EMC 1, 2) = 

EMC 2, 1) = 

EMC 2, 2) = 

TRANSFER MATRIX 

EMC 1, 1) = 

EMC 1, 2) = 

EMC 2, 1) = 

EMC 2, 2) = 

TRANSFER MATRI X 

EMC 1, 1) = 

EMC 1, 2) = 

EMC 2, 1) = 

EMC 2, 2) = 

TRANSFER MATRI X 

EMC 1, 1) = 

EMC 1, 2) = 

EMC 1, 1) = 

EMC 2, 2) = 



X(t) - A X(t) + B XCt - .5) 
+ D x U(t) + D 2 U(t - .5) 



s] 



o 

-5.24 



D l" 



,1)=--,D2(2,1)=--) 



TRANSITION MATRI CES 



00 
00 
00 

00 

00 



PHIC 0) 

. 778801E 

.172270E 

.OOOOOOE 

.606531E 

PHIC 1) 

-.235067E 

- . 1878G6E-01 

-.180539E 01 

-.216281E 00 

PHIC 2) 

. 126127E-01 

.614986E-03 

.196883E 00 

.119 37 7E-01 

PH I ( 3 ) 

-.273338E-03 

-.95 88U8E-0 5 

- . 6Ui+50GE-02 

-.263750E-03 

PHIC k) 

• 31838(+E-05 

.8721U8E-07 

. 100487E-03 

. 309662E-05 

PHIC 5) 

-.231099E-07 

-.519 268E-09 

-.9H011E-06 

-.225906E-07 



TRANSFER MATRIX DELTAC 0) 

) - .256388E 00 



DELC 
DELC 



1, 
2, 



TRANSFER MATRIX DELTAC 1) 

L) = -. 132S18E-01 



DELC 

DELC 



1, 
2, 



TRANSFER MATR 
DELC 1, 



DELC 



2, 



TRANSFER MATR 
DELC 1, 



DELC 



2, 



TRANSFER MATR 
DELC 1, 



DELC 



2, 



TRANSFER MATR 
DELC 1, 

DELC 2, 



) = 



. 20617SE 01 



) = -.210165E 00 

IX DELTAC 2) 

) = .283552E-03 



) = 



.672861E-02 



IX DELTAC 3) 

) = -.327556E-05 

) = -.103763E-05 

IX DELTAC k) 

) = .236515E-07 



) = 



.937b52E-06 



IX DELTAC 5) 

) = -.116723E-09 

) = -.555865E-08 
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TRANSFER MATRIX PH I ( 5) TRANSFER MATRIX DELTA( 6) 

EMC 1, 1) = . lll+£+63E-09 DEL( 1, 1) = . 418U03E- 12 

EM( 1, 2) = .218008E-11 

EM< 2, 1) = . 54U192E-08 DEL( 2, 1) = .232657E-10 

EMC 2, 2) = .112283E-09 

GIVE THE INITIAL STATE (X(l,l)=--- ) 

x(l,l)=O.,0.» 

TIME X(l) X(2) 



.50 


. 2 564E 


00 


.20B2E 


01 


1.00 


. 7980E 


00 


.3102E 


01 


1.50 


.1300E 


01 


.2831E 


01 


2.00 


.1502E 


01 


.1539E 


01 


2.50 


.1332E 


01 


. 2I+06E- 


•01 


3.00 


.9204E 


00 


-.8884E 


00 


3.50 


.5132E 


00 


-.7847E 


00 


4.00 


.3246E 


00 


.1653E 


00 


4.50 


.U295E 


00 


.1364E 


01 


5.00 


. 7391E 


00 


.2150E 


01 


5.50 


. 1067E 


01 


•2155E 


01 


6.00 


. 1237E 


01 


. 1465E 


01 


6. 50 


. 1179E 


01 


.5227E 


00 


7.00 


.9475E 


00 


-.1450E 


00 


7.50 


.6857E 


00 


-.2169E 


00 


8. 00 


.5349E 


00 


.2774E 


00 


8. 50 


.5622E 


00 


.1013E 


01 


9.00 


.7330E 


00 


.1573E 


01 


9. 50 


. 9 4 8 E 


00 


.1683E 


01 


10.00 


.1072F 


01 


.1335E 


01 


10.50 


. 10b5E 


01 


.7646E 


00 


11.00 


.9404E 


00 


•2988E 


00 


11.50 


•7766E 


00 


. 1715E 


00 


12.00 


.6646E 


00 


.4115E 


00 


12.50 


.6579E 


00 


.8506E 


00 


13.00 


.748QE 


00 


. 1234E 


01 


13.50 


.8763E 


00 


.1366E 


01 


14.00 


.9708E 


00 


.1205E 


01 


lit. 50 


.9854E 


00 


.8692E 


00 


15.00 


.9214E 


00 


.5560E 


00 


END OF EXECUTION 








TO CONTINUE, 


GO TO THE TOP 


OF A NEW 


PAGE 


AND PRINT AN 


ASTERI SK 









Figure 5.8 Console transaction for example 5,2,2 
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GIVE ORDER OF SYSTEM (M = ) 

SAMPLING TIME (T = ), FINAL TIME (TF = )| 

m=2, t=.5, tf=ll.* ! 

IS THERE ANY DISTURBING SIGNAL 
yes 

GIVE NUMBER OF INPUT SIGNALS (R = ) 
r = l* 

GIVE THE A MATRIX ( A( 1, 1 ) = --, A( 2 / 1 )=-- ) 
3(1,1)"-. 5, .5,0.* 
3(2,1)=-!. 85,-1. ,1.85* 
a(3,l)=0.,0.,0.* 

GIVE INITIAL STATE (X(l)=--) 
x(l)=0.,0.* 

TERMS OF THE MATRIX EXPONENTIAL 



V(t) - A V(t) 



V(t) - 



X(l>" 
X(2) 
U(l) 



A - 



-.5 .5 

-1.85 -1 





1. 



X(O) 



-[:] 



EM( 


1, 1) = 




697370E 00 




EM( 


1, 2) * 




165714E 00 




EMC 


1, 3) = 




888757E-01 




EM( 


2, 1) = 


-. 


6131U1E 00 




EMC 


2, 2) * 




531656E 00 




EMC 


2, 3) = 




702016E 00 




EMC 


3, 1) = 




000000E 00 




EMC 


3, 2) = 




000000E 00 




EMC 


3, 3) = 


i! 


O0O000E 00 




TIME 


XC1) 




X(2) 




.50 


•888757E- 


■01 


.702016E 


00 


1.00 


.267189E 


00 


.102075E 


01 


1.50 


.i+U<+358E 


00 


.108088E 


01 


2.00 


.57787UE 


00 


.100U22E 


01 


2.50 


.658281E 


00 


.881598E 


00 


3.00 


.69l*033E 


00 


.767104E 


00 


3.50 


.699993E 


00 


.68^312E 


00 


I*. 00 


.690I+29E 


00 


.6366^1E 


00 


(+.50 


.675860E 


00 


.617160E 


00 


5.00 


. 662U72E 


00 


.615736E 


00 


5.50 


.652899E 


00 


.623187E 


00 


6.00 


.6U7^58E 


00 


.633019E 


00 


6.50 


.6U5293E 


00 


.6U1581E 


00 


7.00 


. 6U5202E 


00 


.6U71+61E 


00 


7.50 


,6i*6113E 


00 


.6506U3E 


00 


8.00 


.61+7276E 


00 


.651776E 


00 


8.50 


.6U827UE 


00 


.651666E 


00 


9.00 


.6I+8953E 


00 


.650995E 


00 


9.50 


.6^931UE 


00 


.650222E 


00 


10.00 


.6i|9i+38E 


00 


.6U9590E 


00 


10.50 


.649420E 


00 


.6**9178E 


00 


11.00 


.649339E 


00 


.61+8970E 


00 


END OF EXECUTION 








TO CONTINUE, 


GO TO THE TOP OF A NEW PAGE 




AND PRINT AN 


ASTERISK 









Figure 5,10 Console transaction for example 5.2.3 
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GIVE ORDER OF SYSTEM (M = ) 

DESIRED SAMPLING TIME (T = ) 

TIME DELAY C TO = ), FINAL TIME (TF = ) 

m=2,t=.5,td=.5, tf=15.* 

IS THERE ANY DISTURBING SIGNAL 
yes 

NUMBER OF INPUT SIGNALS (R = ) 
r = l* 

GIVE THE A MATRIX < A( 1, 1 ) =--, A( 2, 1 ) =-- ) 

a(l, l)=-.5,.5* 

a(2,l)»0.,-l.* 

GIVE THE B MATRIX ( B ( 1, 1 ) = --, B ( 2, 1 ) =-- ) 

b(l,l)=0.,0.* 

b(2,l) — 1.85,0.* 

GIVE THE Dl MATRIX ( Dl( 1, 1) =--, Dl( 2, 1 ) =-- ) 

dl(l,l)=0.* 

dlC 2, 1) =1.85* 

GIVE THE D2 MATRIX ( D2( 1, 1 ) =--, D2 ( 2, 1 ) =-- ) 

d2(l,l)=0.* 

d2(2,l)=0.* 

DO YOU WISH TO HAVE THE TRANSITION MATRICES 

no 

GIVE THE INITIAL STATE (X(l,l)= ) 

x(l,l)=0.,0.* 

TIME X(l) X(2) 
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X(t) - A X(t) + B X(t - .5) 
+ DjUCt) + D 2 U(t - .5) 



[,3S] 

[o] 



D l- 



D 2- 



.50 
1.00 

1.50 



00 
50 
00 

50 
00 
50 
00 



5.50 

6.00 

6.50 

7.00 

7.50 

8.00 

8.50 

9.00 

9.50 

10.00 

10.50 

11.00 

11.50 

12.00 

12. 50 
13.00 

13. 50 
14.00 
14.50 
15.00 



.9052E-01 

.2848E 00 

.4952E 00 

.6644E 00 

.7664E 00 

.8015E 00 

.7862E 00 

.7431E 00 

.6932E 00 

.6512E 00 

.6245E 00 

.6138E 00 

.6158E 00 

.6251E 00 

.6369E 00 

•6U72E 00 

.6541E 00 

.6572E 00 

.6572E 00 

.6552E 00 

.6525E 00 

.6499E 00 

.6U82E 00 

.6473E 00 

.6U72E 00 

.6476E 00 

.6482E 00 

.6488E 00 

.6493E 00 

.6495E 00 



.7279E 00 
.1U3E 01 
.1282E 01 
•121UE 01 
.103UE 01 
.8266E 00 
.6539E 00 
.5448E 00 
.5021E 00 
.5114E 00 
.5508E 00 
•5995E 00 
.6I+21E 00 
.670UE 00 
.6829E 00 
.6825E 00 
.6740E 00 
.6627E 00 
.6523E 00 
.6U50E 00 
.641UE 00 
.6411E 00 
.6429E 00 
.6U55E 00 
.6U80E 00 
.6U99E 00 
.6508E 00 
.6510E 00 
.6507E 00 
.6501E 00 



Figure 5.12 

Console transaction for example 
5.2.4 
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p-t <N ro 
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loadgo timdel delfor pertur 

W 12 50.1 

EXECUTION. 

GIVE ORDER OF SYSTEM (M = ) 

DESIRED SAMPLING TIME (T = ) 

TIME DELAY CTD = ), FINAL TIME (TF = } 

m=2, t=.5, td=l., tf =15.* 

IS THERE ANY DISTURBING SIGNAL 
yes 

NUMBER OF INPUT SIGNALS (R = ) 
r~l* 

GIVE THE A MATRIX ( A C 1, 1 ) =--, A( 2, 1 ) =-- ) 

3(1,1)=-. 5, .5* 
a(2,l)=0.,-l.* 

GIVE THE B MATRIX { B{ 1, 1) =--, B{ 2, 1) =-- ) 

b(l,l)=0.,0.* 

b(2, 1)=-1.85 # 0. * 

GIVE THE Dl MATRIX (Dl< 1, 1) =--, Dl< 2, 1 ) =-- ) 

dl(l / l)=0.* 

dlC2,l)=1.85* 



GIVE THE 02 MATRIX (D2 

d2(l,l)=0.* 

d2(2, 1)=0.* 

DO YOU WISH TO HAVE TH 

yes 

TRANSFER MATR 

EMC 1, 1) 

EMC 1, 2) 

EMC 2, 1) 

EMC 2 , 2 ) 

TRANSFER MATR 

EMC 1, 1) 

EMC 1, 2) 

EMC 2, 1) 

EMC 2, 2) 

TRANSFER MATR 

EMC 1, 1) 

EMC 1, 2) 

EM ( 2 , 1 ) 

EMC 2, 2) 

TRANSFER MATR 

EMC 1, 1) 

EMC 1, 2) 

EMC 2, 1) 

EMC 2 , 2 ) 

TRANSFER MATR 

EMC 1, 1) 

EMC 1, 2) 

EMC 2, 1) 

EMC 2, 2) 

TRANSFER MATR 

EMC 1, 1) 

EMC 1, 2) 

EMC 2 , 1 ) 

EMC 2, 2) 



X(t) - A X(t) + B X(t - 1) 
+ D x U(t) + D 2 U(t - 1) 






(l,l)=--,D2(2,l)=--> 

E TRANSITION MATRICES 

IX PH I ( ) 

.773801E 00 

. 172270E 00 

.OO0O00E 00 

.606531E 00 

IX PHIC 1) 

= -.829913E-01 

= -.663267E-02 

= -.637399E 00 

= -.763586E-Q1 

IX PHIC 2) 

.157213E-02 

.766560E-0I* 

. 2f*5t*09E-01 

. 1U95U8E-02 

IX PH I ( 3 ) 

= -.120288E-0U 

= -.1+21960E-06 

= -.283627E-03 

= -.116068E-OI* 

IX PH I ( U ) 

.1+94666E-07 

.135504E-08 

.156125E-05 

.481116E-07 

IX PH I C 5 ) 

= -.126765E-09 

= -.284835E-11 

= -.501365E-08 

= -.123917E-09 






TRANSFER MATRIX DELTAC 0) 

L) = .905188E-01 



DELC 


1 




DELC 


2 




TRANSFER 
DELC 1 


MATR 


DELC 


2 




TRANSFER 
DELC 1, 


MATR 


DELC 


2, 




TRANSFER 
DELC 1, 


MATR 


DELC 


2, 




TRANSFER 
DELC 1, 


MATR 


DELC 


2, 




TRANSFER 
DELC 1, 


MATR 


DELC 


2, 





) = 



•727918E 00 



IX DELTAC 1) 

) = -.165553E-02 

) = -. 261961+E-01 

IX DELTAC 2) 

) = .12i*783E-0U 



) = 



.296106E-03 



IX DELTAC 3) 

) = -.508918E-07 

) = -.16121«»E-05 

IX DELTAC it) 

) = .129736E-09 



) = 



.5U338E-08 



IX DELTAC 5) 

) = -. 226048E-12 

) = -. 1076U9E-10 
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GIVE THE INITIAL STATE (X(l,l)=---) 
x(l,l)=0.,0.* 

TIME X(l) X(2) 



.50 


.9052E-01 


.7279E 00 


1.00 


.2864E 00 


.1169E 01 


1.50 


.5134E 00 


.1411E 01 


2.00 


. 7191+E 00 


. 1UU4E 01 


2.50 


.8664E 00 


.1306E 01 


3.00 


.9369E 00 


.1063E 01 


3.50 


.9328E 00 


.7864E 00 


i+.OO 


.8712E 00 


. 5U17E 00 


4.50 


.7771E 00 


.3720E 00 


5.00 


.6770E 00 


.2960E 00 


5.50 


.5328E 00 


.3093E 00 


6.00 


. 538UE 00 


.3897E 00 


6.50 


.5185E 00 


.5062E 00 


7.00 


.5299E 00 


.6269E 00 


7.50 


.5634E 00 


.7262E 00 


S.00 


.6073E 00 


.7881E 00 


8.50 


.6503E 00 


.8080E 00 


9.00 


.6838E 00 


.7909E 00 


9.50 


.7029E 00 


•7482E 00 


10.00 


.7068E 00 


.6944E 00 


10.50 


.6980E 00 


. 6U29E 00 


11.00 


.6810E 00 


.6038E 00 


11.50 


.6611E 00 


.5825E 00 


12.00 


.6430E 00 


.5795E 00 


12.50 


.6301E 00 


.5915E 00 


13.00 


.6239E 00 


.6128E 00 


13.50 


.6243E 00 


.6370E 00 


14.00 


.6296E 00 


.6584E 00 


14.50 


.6379E 00 


.6733E 00 


15.00 


.6466E 00 


.6800E 00 


END OF EXECUTION 




TO CONTINUE, 


GO TO THE TOP 


OF A NEW PAGE 


AND PRINT AN 


ASTERISK 





Figure 5.14 Console transaction for example 5.2.5 
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X_Ct) 



"o 


X 


X(t) + 


,2 





XCt - T) + 





-I 










-.1 




1 



U(t - T) (5.26) 



It Is assumed the sampling time equal to the time delay* That is, 



T - T - 



■ min 



(5.27) 



Koepcke reported the following results of the plant transition 
matrices and the control transition matrices i 



*o 



*! - 



*2 - 



*3 " 



** - 



* 5 - 



.7071068 .7071068 

-.7071068 .7071068 

" .1338340 .0277680" 

-.0277680 -.0782980 

" .0109582 .0022524 

-.0022524 .0026278 

~ .0005903 .0000742" 

-.0000742 -.0000854 

" .0000236 .0000026 

-.0000026 .0000013 

~ .0000008 .0000001 

-.0000001 -.0000000 



A 3 - 



A 5 - 



.00000 
.00000 

* .2928932" 
.7071068 

" .0075873" 
-.0308106 

~ .0004532" 
.0007349 

' .0000119" 
_-. 0000165, 

" .0000003 
.0000002 



The time response of the system was obtained assuming a step 
input and zero initial conditions for the integrators. 

The solution is depicted In figures 5.18 and 5.19. 
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In a similar way, this same example was tested assuming no lags 
in the system, that is 



V(t) 



.2 1 
-1 -.1 




V(t) 



(5.28) 



where 



V(t) - 



x 2 (t) 
u(t) 



(5.29) 



The evaluation of the state is shown in figures 5.16 and 5.17. 

It is interesting to compare the transient response in both cases. 
As it can be seen in the plots (figures 5.17 and 5.19), the case with 
delay is something less unstable than the linear one with delay equal to 
zero. 



68 



GIVE ORDER OF SYSTEM CM = ) 

SAMPLING TIME (T = ), FINAL TIME (TF = ) ! 

m=2, t=.5, tf=15.* | 



IS THERE ANY DISTURBING SIGNAL 
yes 

GIVE NUMBER OF INPUT SIGNALS (R = ) 
r = l* 



} V(t) - A V(t) 



! V(t) - 



GIVE THE A MATRIX ( AC 1, 1 ) =--, AC 2, 1 ) =— ) ! 

aCl,l)=.2,l.,0.* i 

a(2,l)=-l.,-.l,l.* |A 

a(3,l)=0.,0.,0.* j 



X(l) 
X(2) 
U(l) 



.2 


1 





-1 


-.1 


1 












GIVE INITIAL 


STATE CXC1) 




» i 




xCl)=0.,0.* 






1 

! 

1 




TERMS OF THE 


MATRIX EXPONENTIAL 




EMC 


1, 1) = 




976370E 00 




EMC 


1, 2) * 




U92031E 00 




EMC 


1, 3) = 




121+527E 00 




EMC 


2, 1) = 


-, 


I+92031E 00 




EMC 


2, 2) = 




828760E 00 




EMC 


2, 3) ■ 




1+67126E 00 




EMC 


3, 1) = 




000000E 00 




EMC 


3, 2) = 




000000E 00 




EMC 


3, 3) = 


i! 


OOOOOOE 00 




TIME 


XC1) 




XC2) 




.50 


. 12U527E 


00 


.U67126E 


00 


1.00 


.1+75952E 


00 


. 792990E 


00 


1.50 


.979t*07E 


00 


.8901U1E 


00 


2.00 


.151877E 


01 


. 7229<+0E 


00 


2.50 


.196311E 


01 


. 318989E 


00 


3.00 


.219820E 


01 


- . 23i*422E 


00 


3.50 


. 215 5U1+E 


01 


-.808739E 


00 


1*. 00 


.183111E 


01 


-.126367E 


01 


k.SO 


.129060E 


01 


-.1U8112E 


01 


5.00 


.655877E 


00 


-.139538E 


01 


5.50 


. 783336E- 


-01 


-.101203E 


01 


6.00 


-.296938E 


00 


-. U101U3E 


00 


6.50 


-. 367197E 


00 


. 273318E 


00 


7.00 


-.995121+E- 


-01 


.87U313E 


00 


7.50 


.d57555E 


00 


.124068E 


01 


8.00 


. 118173E 


01 


. 127022E 


01 


8.50 


.190332E 


01 


.938392E 


00 


9.00 


.24I+U59E 


01 


.308336E 


00 


9.50 


.266306E 


01 


-.1+80150E 


00 


10.00 


,2i*88UE 


01 


-.12U111E 


01 


10.50 


.19U3U7E 


01 


-.178583E 


01 


11.00 


.11U338E 


01 


-.196915E 


01 


11.50 


.272011E 


00 


-.1727U0E 


01 


12.00 


-.U59826E 


00 


-. 109832E 


01 



X(O) 



■m 
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To 
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jie 5.16 Console transaction ior example 3.3 
when time delay « U 
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GIVE ORDER OF SYSTEM (M = ) 

DESIRED SAMPLING TIME (T = ) 

TIME DELAY CTD = ), FINAL TIME CTF = ) 

in=2, t = . 78 5 39 8 2, td = . 78 5398 2, tf = 2 0.* 

IS THERE ANY DISTURBING SIGNAL 
yes 

NUMBER OF INPUT SIGNALS CR = ) 
r = l* 

GIVE THE A MATRIX C AC 1, 1 ) = --, AC 2, 1 ) =-- ) 

a(l,l)=0., 1.* 
aC2,l)=-l.,0.* 

GIVE THE B MATRIX ( B ( 1, 1 ) =--, B ( 2, 1 ) =-- ) 

b(l,l) = .2,0.* 

bC2,l)=0.,-.l* 

GIVE THE Dl MATRIX (D 1C 1, 1 ) =--, Dl C 2, 1 ) = -- ) 

dl(l,l)=Q.* 

dlC 2, 1)=0.* 



X(t) - A X(t) + B X(t - j) 



+ D x U(t) + D 2 U(t - ~) 



A - 





-1 



J] 



[i -',] 



GIVE THE D2 MATR 1 X CD2C 1, 


1)=--,D2.<2,1)=- 


-) j 


d2(l, 1)=0.* 
d2(2, 1)=1.* 






• ^ [J] 

CES 


DO YOU WISH TO HAVE THE TRANSITION MATRI 


yes 








TRANSFER 


MATRIX 


PH 1 C ) 


TRANSFER MATRIX DELTAC 0) 


EMC 1, 


1) = 


.707107E 00 


DELC 1, 1) = .00000OE 00 


EMC 1, 


2) = 


. 707107E 00 




EMC 2 , 


1) = 


-.707107E 00 


DELC 2, 1) - .OOOOOOE 00 


EMC 2, 


2) = 


.707107E 00 




TRANSFER 


MATRIX 


PH 1 ( 1 ) 


TRANSFER MATRIX DELTAC 1) 


EMC 1, 


1) = 


. 13383!jE 00 


DELC 1, 1) = .292893E 00 


EMC l, 


2) = 


.2776K0E-01 




EM C 2 , 


1) = 


-.277680E-01 


DELC 2, 1) = . 707107E 00 


EMC 2, 


2) = 


-.782980E-01 




TRANSFER 


MATRIX 


PHI ( 2) 


TRANSFER MATRIX DELTAC 2) 


EMC 1, 


1) = 


.109582E-01 


DELC 1, 1) = .758732E-02 


EMC 1, 


2) = 


.225237E-02 




EMC 2, 


1) = 


-.225237E-02 


DELC 2, 1) = -.308106E-01 


EMC 2, 


2) = 


. 262 783E-02 




TRANSFER 


MATRIX 


PHIC 3) 


TRANSFER MATRIX DELTAC 3) 


EMC 1, 


1) = 


.5903U3E-03 


DELC 1, 1) = .[+53237E-03 


EMC 1, 


2) = 


. 7U1 76 5E-04 




EMC 2, 


1) = 


-, 7<U765E-0i4 


DELC 2, 1) = . 73U9 07E-03 


EMC 2, 


2) = 


-.853680E-0if 




TRANSFER 


MATRI X 


PHIC k) 


TRANSFER MATRIX DELTAC k) 


EMC 1, 


1) = 


.235559E-0U 


DELC 1, 1) = .U8773E-0I* 


EMC 1, 


2) = 


.25925UE-05 




EMC 2, 


1) = 


- . 2 5925i+E-05 


DELC 2, 1) = - . 16U 710E-0I* 


EMC 2, 


2) = 


. 130299E-05 




TRANSFER 


MATRI X 


PHIC 5) 


TRANSFER MATRIX DELTAC 5) 


EMC 1, 


1) = 


. 71+3663E-Q6 


DELC 1, 1) = .3U1U10E-06 


EMC 1, 


2) = 


.65U65E-07 




EMC 2, 


1) = 


-.65U65E-07 


DELC 2, 1) = .217079E-06 


EMC 2, 


2) = 


-.290989E-07 
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TRANSFER MATRIX PH I ( 6) TRANSFER MATRIX DELTA( 6) 

EMC 1, 1) = . 19762i*E-07 DEL( 1, 1) = .739099E-08 

EM( 1, 2) = .150533E-08 

EMC 2, 1) = -.150533E-08 DELC 2, 1) = -.367553E-08 

EMC 2, 2) = .218458E-09 

GIVE THE INITIAL STATE ( XC 1, 1 ) = ) 

x(l, 1>=0.,0.* 

TIME XC1) X(2) 



.79 


.OOOOE 00 


.OOOOE 00 


1.57 


.2929E 00 


.7071E 00 


2,36 


.1008E 01 


.9692F. 00 


3.11* 


. 175SE 01 


•5864E 00 


3.93 


•2125E 01 


-.2538E 00 


d.71 


.1889E 01 


-.1100E 01 


5.50 


.1158E 01 


-.1478E 01 


6.28 


•3207E 00 


-.1159E 01 


7.07 


-.1582E 00 


-.2928E 00 


7.85 


.3256E-02 


.6571E 00 


8.64 


.7401E 00 


.1163E 01 


9.1+2 


.1663E 01 


.9241F 00 


10.21 


. 2263R 01 


. UU67E-Q1 


11.00 


.2192E 01 


-.1009E 01 


11.78 


. H+62E 01 


-.1654E 01 


12,57 


.4567E 00 


-.1514E 01 


13.35 


-.2735E 00 


-.6351E 00 


11*. 11+ 


-.3088E 00 


,5191+E 00 


14.92 


.39 80E 00 


. 1315E 01 


15. 71 


.1481E 01 


.1292E 01 


16.1*9 


.2349E 01 


,t*317E 00 


17.28 


. 2 5 9 E 01 


-.8138E 00 


18.06 


.1842E 01 


-.1775E 01 


18.85 


•6889E 00 


-.1890E 01 


19.63 


-.3244E 00 


-. 1067E 01 


20.42 


-.6256E 00 


.2718E 00 


END OF EXECUTION 




TO CONTINUE, 


GO TU THE TOP 


OF A NEW PAGE 


AND PRINT AN 


ASTERISK 





Figure 5.18 Console transaction for example 5,3 
when time delay - — min 
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CHAPTER 6 
COMMENTS AND SUGGESTIONS FOR FUTURE RESEARCH 

AT 
In obtaining e by the use of a digital computer the virtues of 

the series expansion technique are its simplicity and ease in programming. 

It is not necessary to find the eigenvalues of A. There is, however, some 

computational disadvantage to the series expansion method. This comes from 

the convergence requirements for the series. In general, it is reasonable 

AT 
to compute e by the power series when T is small. The running time for 

the matrix exponential simulation will be among the longest of various 

schemes. Use of the Jordan Canonical form, for example, requires 

considerably more programming, but will run in a fraction of time needed 

for the series solution. 

Some suggestions concerning the bound on the error In the 
evaluation of the matrix exponential when the matrix A is known with some 
error are given by Levis (10). 

The simulation technique for linear time- invariant dynamic 
systems has been tested, and it was found that the use of the augmented A 
matrix (X(t) - A X(t) + D U(t) can be expressed as V(t) - £ ° V(t) , where 
V(t) - \~\ ) greatly improved the procedure. The reason is that the actual 
reduction of the elements of the augmented matrix times T to values less 
than one can be performed successfully. However, this method cannot be 
used for calculating the digital version of the control transition matrix. 

Another scheme that can be used to check the error bound in the 
state is to divide the time region of interest in two or three parts. 
Preferably these times should be powers of two times the sampling time. 
Next, compute the matrix exponential at the desired sampling time. 



75 



Recursively multiply it until the matrix exponential is found for the other 
selected times. The state at those times can be found and saved. Nov, 
using the recursive process of state evaluation at the sampling time, 
compare the state with the selected ones. If the error is unacceptable, 
the state with less error can be used as a new initial condition, and the 
procedure may be continued. 

It was found in chapter 3 that the elements C^ . form an array 
of Infinite order. The first row is of main Importance because its elements 
are the terms of e AT . Therefore, the truncation technique already discussed 
can be used. 

In a similar fashion, the elements C^ ^ are actually the terms of 
the infinite series e BT . It is reasonable to expect smaller values of these 
norms as "i" grows. Therefore, Intuitively the number of terms used to 
truncate the first row can be used to truncate $i(t), *2(t), etc. It would 
be interesting to make a study about how the truncation terms should be 
taken in each row in order to save computation time while maintaining 
accuracy. 
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TRANS 



Purpose: to compute the time response of linear time- invariant 
By stems. 

Inputs: order of system (M - ) ; sampling time (T - ) ; final 
time (TF - ) ; number of input signals (R - ) ; the 
augmented A matrix and the initial state (X(l) - ). 

Outputs: the transition matrix; the current time; and the state 
of the system. 

Remarks: main program. Subroutine called by TRANS: EXPMAT, and 
D1STUR. 
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PROGRAM COMMON A, EM, M, RIJ, R, x 

DIMFNSION X(2»)»Y(2u),E(2U)»PF(20)tXK2 0) 

DIMENSION FMPI40CH) ,A< 40'J,H ) »EM ( 400,H) 

INTEGER I ,J,M,R,wiSH 

FORMAT VARIAflLF FM 

VECTOR VALUES H=2»1»Q 
MAGDA PRINT COMMENT SGIVE ORDER OF SYSTEM (M - ) S 

PRINT COMMENT iSAMPLING TIME (T = ), FINAL TIME ( TF = )S 

READ DATA 

PRINT COMMENT i $ 

FM = M 

PRINT COMMENT SIS THERE ANY DISTURBING SIGNAL* 

READ FORMAT S3 .WISH 

VECTOR VALUFS S3 = $ C3*$ 

WHENEVER WI5H.F.SYES$ 

PRINT COMMENT $> i 

PRINT COMMENT £GIVE NUMBER OF INPUT SIGNALS (R = )$ 

READ DATA 

M = M+R 

OTHERWISE 

R = 

END OF CONDITIONAL 
H(?]=M 

PRINT COMMENT S £ 

PRINT COMMENT SGIVE THE A MATRIX (AU,1> = — ,A(2,1)= — )$ 
THROUGH LUPE, FOR 1 = 1,1, I. G.M 
LUPE READ DATA 

PRINT COMMENT S 3 

PRINT COMMENT SGIVE INITIAL STATE (X(l)=— IS 
READ DATA 

THROUGH ALICIA, FOR I = 1 » 1 . I . G. ( M-R ) 
ALICIA XI ( I )=X( I) 
TA = T 

WHENEVFR R.NF.U 
EXECUTE DISTUR. (TA) 
j=M-R+i 

THROUGH JULIA, FOR I=J,1,I.G.M 
XI ( J )=X( J) 
JULIA CONTINUF 

END OF CONDITIONAL 
THROUGH ALMA, FOR I - 1 , 1 , I . G. ( M ) 
ALMA E(I)=0. 
TZ = T 

EXECUTE ^XPMAT.CT) 
THROUGH FANNY, FOR I=l»liI.G.M 
THROUGH FANNY, FOR J=1,1,J.G.M 
PRINT FORMAT CUA TRO , I , J , EM ( I , J ) 

VECTOR VALUES CUATRO = S 1H , S8 , 3HEM ( , I 4 , 1H , , I A , 3H ) =»E14.6*1 
FANNY CONTINUF 

THROUGH MARTA, FOR I = 1 , 1 , I .G . ( M-R ) 
THROUGH MARTA, FOR J=1,1,J.G.M 
MARTA EMP( I ,J) = EM( I ,J I 

WHENEVER (M-RJ.L.6 
PRINT COMMFNT $ & 

PRINT FORMAT SI, ( I = 1 , 1 , I .G. < M-R ) , I ) 

VECTOR VALUES SI = i » S6 , 4HT I ME , SB , ■ FM • ( 2HX ( , I 1 , 1H ) , S 1 2 ) /*$ 
END OF CONDITIONAL 
TRANSFER TO TERESA 
OLGA TA=TA+TZ 
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TFRFSA 



MAR I A 
FLFNA 



ROSA 



F.STHFR 
ROSANA 



SARA 
CARMEN 

I. ILIA 



WHFNFVFR R.MF.O 

EXECUTF DISTUR. (TA) 

FND OF CONDITIONAL 

THROUGH ELENA, FOR I = 1 » 1 » I .G. ( M-R > 

PE( I )=0. 

Y( I )=0. 

THROUGH MARIA* FOR J=1»1»J.6.M 

Y{ I ) =Y( I )+FMP< I ,J)*X< J ) 

PF( I )= ( pM-P{ I,J)+RIJ)*t:tJ)+RIJ*X(J)+PE( I ) 

CONT INUF 

CONT INUF 

f nj o R m = n . 

THROUGH ROSA, FOR I = 1 » 1 » I . G. ( M-R ) 

ENOR M=ENORM+. ABS . ( PF{ I ) ) 

WHENEVER <FNORM.GF. ( 10. .P. -07) 

T = TA 

EXFCUTE EXP^AT.(T) 

THROUGH ROSANA, FOR I = 1 » 1 » I .G . ( M-R ) 

P E ( I ) = . 

Y( I ) =0. 

THROUGH ESTHER. FOR J=1,1,J.G.M 

Y { I ) = Y ( I )+EM( I »J)*XI ( J ) 

PE( I ]=PF( I ]+RI J*XI < J) 

CONT INU r 

CONTINUF 

OTHERWISE 

TRANSFFR TO SARA 

FND OF CONDITIONAL 

THROUGH CARMEN, FOR I = 1 , 1 , I . G. ( M-R ) 

X( I )=Y( I ) 

F( I )=PF E I ) 

j=M-R+l 

THROUGH LILIA, FOR I=J,1,I.G.M 

F( J)=0. 

WHFNFVFR ( m-R ) .L.6 

PRINT FORMAT 52 » T A , X ( 1 ) . . . X ( M-R > 

VECTOR VALUES S2 = i » 54 » F6 . 2 » ' FM " ( S3 , E 1 A. 6 ) *i 

OTHERWISE 

PRINT RFSULTS TA 

PRINT RESULTS X ( 1 ) . . . X ( M-R } 

FND OF CONDITIONAL 

WHFNFVFR TA.L.TF, TRANSFFR TO OLGA 

PRINT COMMENT SEND OF EXECUTIONS 

PRINT COMMENT STO CONTINUE, GO TO THE TOP OF A NEW PAGES 

PRINT COMMENT SAND PRINT AN ASTERISKS 

READ DATA 

TRANSFFR TO MAGDA 

END OF PROGRAM 
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EXP MAT 



irpose ; ; ■ compute tiae matrix exponential, 



Remarks' ^ coroutine called by TRANS, 
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EXTERNAL FUNCTION IT) 

PROGRAM COMMON A, EM, M, RIJ, R, X 

DIMENSION Al 400tH) »EM( 4 0U,H> »TFRM<40 0»H) »NTERM< ^0 0»H) 

DIMENSION X< 20) ,B(4uO,H) 

VECTOR VALUES H=2,l,0 

INTEGER K»I,J,L»M,LL,Y»Q 

ENTRY TO EXPMAT. 

H(2) = y 

THROUGH ELENA, FOR 1 = 1,1, I. G.M 

THROUGH ELENA, FOR J=1,1,J.G.M 

B( I »J)=A{ I .J) 
ELENA B( I * J)=R( I »J)*T 

AMIN = R( 1 ,1 ) 

THROUGH DIANA, FOR I=2,l,J.G.v 

WHENEVER P( I , I ) .L.AMIN, AVIN=B(I»I> 
DIANA CONTINUE 

FAC=EXP. ( AMIN) 

THROUGH OLGA* FOR 1 = 1,1, I. G.M 
OLGA B( I ♦ I ) = B< I »I 1-AMIN 

Y=.ABS.B< 1,1) 

THROUGH SARA, FOR 1=1, 1,1. G.M 

THROUGH SARA* FOR U=1»1»U.G.M 

WHENEVER . ABS . ( B ( I , U > ) ,G . ( Y + . ) , Y = . A8S . ( B ( I , J ) ) 
SARA CONTINUF 

TAP=1. 

YE=Y+0. 

THROUGH ALMA, FOR Q=l,l,Q.G.lO 

TAP=2.*TAP 

WHENEVER TAP.GE.YE, TRANSFER TO ESTHER 
ALMA CONTINUF 
ESTHER Y=TAP 

THROUGH YOLIS, FOR 1 = 1, 1,1. G.M 

THROUGH YOLIS, FOR J=1,1,J.G.M 

B( I » J)=B( I ,J) /(Y+0. ) 
YOLIS TFRM( I ,J)=B( I ,J) 

LL = 
GLORIA MAXH=0. 
MAXV=0. 

THROUGH MARIA, FOR 1 = 1, 1,1. G.M 
SUMH=0. 
SUMV=0, 

THROUGH ROSANA, FOR J=1»1,J.G.M 
SUMH=SUMH+.ABS.TERM( I , J) 
SUMV=SUMV+.ABS.TERM< J, I ) 
ROSANA CONTINUE 

WHENEVER SUMH .G.MAXH »MAXH=SUMH 
WHENEVER SUMV.G. MAXV ,MAXV=SUMV 
MARIA CONTINUE 
NORM=MAXH 

WHENEVER MAX V. L . NORM ,NORM = MAXV 
WHENEVER LL.NE.O, TRANSFER TO DELIA 
SOLO=NORM 
K=2.*N0RM 

WHENEVER K.L.2, K=2 
IN=</2 

VECTOR VALUES CINCO = S1H ,2HK=,I4*i 
THROUGH SUSANA, FOR 1 = 1, 1,1. G.M 
THROUGH SUSANA, FOR J=1,1,U.G.M 
UNIT=0. 
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WHENEVER J.E. I » UNIT = 1. 

EMf I ,J)=UNIT+B< I »J> 
SUSANA / CONTINUE 
ISABEL WHENEVER LL.GE.K, TRANSFER TO GLORIA 

LL=LL+1 

THROUGH LILIAt FOR L=1»1»L.G.M 

THROUGH LILIAt FOR 1 = 1,1,1. G.M 

NTERM(L, I )=0. 

THROUGH EVA, FOR J=1»1»J.G.M 

NTERM(L, I )=NTERM(L»I )+B(L»J)*TERM( Jtl ) 
EVA CONTINUE 

EM<L»I)=EM(L»I )+NTERM(L, I ) /(LL+1. ) 
I.ILIA CONTINUE 

THROUGH AURORA, FOR 1-1 • 1» I .G.M 

THROUGH AURORA, FOR J=1»1»J.G.M 
AURORA TERM( I,J>=NTERM(I,J)/(LL+1.) 

TRANSFER TO ISABEL 
DELIA EPS=SOLO/(K+2.) 

RIJ=NORM*SOLO/( <K+1.)*< l.-EPSJ ) 

THROUGH JULIA, FOR 1=1,1, 1 .G.M 

THROUGH JULIA, FOR J=1»1»J.G.M 

WW=.ABS. (EMf I ,J)*lO..P.-7) 

WHENEVER RU.G.WW 

K=K+IN 

TRANSFER TO ISABEL 

OTHERWISE 

TRANSFER TO JULIA 

END OF CONDITIONAL 
JULIA CONTINUE 

THROUGH ALICIA, FOR LL= 1 , 1 , LL .G.Q 

THROUGH MARTA, FOR L=1»1»L.G.M 

THROUGH MARTA, FOR 1 = 1,1,1. G.M 

TERM(L,I )-0. 

THROUGH MAGUE, FOR J=1»1»J.G.M 

TERM (L, I ) -TERM (L.I ) +EM ( L , J ) *EMf J , I ) 
MAGUE CONTINUE 
MARTA CONTINUE 

THROUGH OLIVIA, FOR 1 = 1,1, I. G.M 

THROUGH OLIVIA, FOR J=1»1,J.G.M 
OLIVIA EMf I ,J)=TERM( I » J ) 
ALICIA CONTINUE 

PRINT COMMENT $ $ 

PRINT COMMENT $ TERMS OF THE MATRIX EXPONENTIALS 

THROUGH CARMEN, FOR 1 = 1,1, I. G.M 

THROUGH CARMEN, FOR J=1,1,J.G.M 
CARMEN EM( I,J)=FAC*EM< I.J) 

FUNCTION RETURN 

END OF FUNCTION 



at the current 
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TIMDEL 

Purpose: to compute the time response of linear systems with 
lumped parameters and time delays. 

Inputs: order of system (M - ) ; sampling time (T » ) ; time 
delay (TD - ) ; final time (TF - ) ; number of input 
signals (R - ) ; the A matrix; the B matrix; the D 
matrix; the D matrix; the initial state (X(l,l) - ). 

Outputs: the plant transition matrices, the control transition 
matrices if desired; the current time; and the state 
of the system* 

Remarks: main program. Subroutines called by TIMDEL: DELFOR, 
and PERTUR. 
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PROGRAM COMMON EM , DELF »M »R »W * A , B »D1 > D2 »U 

DIMENSION FM(40O0,H},nFLF(4GGO,H),X(4OO»G).A(4OQ»G> 

DIMENSION B ( 400, G) »D1( 400, G) . D2( 400, G) »U< 40 0,E> 

INTEGER I,J,K,L,M,N,LL,Z,k*,',R,REL»MM,WlSH,JJ 

FORMAT VARIABLE FM 

VECTOR VALUFS G=2»l»0 

VECTOR VALUES E=2,l»0 

VECTOR VALUES H=3. 1.0.0 
MAGDA PRINT COMMENT SGIVE ORDER OF SYSTEM (M = >$ 

PRINT COMMENT SDESIRED SAMPLING TIME (T = )$ 

PRINT COMMENT JTIME DELAY ( TD = >t FINAL TIME ( TF = ) S> 

RFAD DATA 

FM = M 

PRINT COMMENT 5 $ 

PRINT COMMENT SIS THERE ANY DISTURBING SIGNALS 

READ FORMAT S3.WISH 

VECTOR VALUES S3 - S C3*£ 

WHENEVER WISH.E.JYESS 

PRINT COMMENT S S 

PRINT COMMENT SNUMBER OF INPUT SIGNALS (R = )$ 

RFAD DATA 

OTHERWISE 

R = 

END OF CONDITIONAL 

REL=TD/T+0.2 

G( 2 ) =M 

H( ? ) =M 

H(3)=M 

PRINT COMMENT $ S 

PRINT COMMENT SGIVE THE A MATRIX (A(l»l)= — »A<2»1)= — )$ 

THROUGH MELA, FOR 1 = 1.1.1. G.M 
MELA READ DATA 

PRINT COMMENT £ S 

PRINT COMMENT SGIVE THE B MATRIX <B<1»1>= — ,R<2»1>=— )S 

THROUGH MALFNA, FOR 1=1.1, 2. G.M 
MALFNA RFAD DATA 

WHENEVER R . E. , TRANSFER TO JULIA 

PRINT COMMFNT S S 

PRINT COMMENT SGIVE THE Dl MATRIX (01(1.1)= — ,DK2.1> = — )* 

THROUGH ALMA. FOR 1=1,1.1. G.M 
ALMA READ DATA 

PRINT COMMENT $ S 

PRINT COMMENT SGIVE THE D? MATRIX (D2<1.1>= — »D2< 2 1 1 > =-- > * 

THROUGH BERTA, FOR 1 = 1.1. I. G.M 
BFRTA RFAD DATA 
JULIA EXFCUTE DELFOR. <T) 

PRINT COMMFNT SDO YOU WISH TO HAVE THE TRANSITION MATRICES* 

READ FORMAT S3. WISH 

WHENEVER WISH.F.SYESS 

THROUGH DULCE, FOR L=1,1,L.G.W 

LL=L-1 

WHENEVER R.E.O 

PRINT FORMAT OCHO.LL 

VECTOR VALUES OCHO = $1H , S8 . 1 5HTRANSFER MATRIX, S2, 
14HPHI ( ,I4,1H)*$ 

OTHERWISE 

PRINT FORMAT SEIS.LL.LL 

VECTOR VALUE5 SEIS= S1H ,58 , 1 5HTR ANSFER MATRIX, S2» 
14HPHI ( .I4.1H) .S8.22HTRANSFER MATRIX DELT A( , I 4 , 1H ) *S 
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CFLIA CONTINUF 

WHENEVER Z.E.I. L=L-1 
GLORIA THROUGH ALICIA, FOR K= 1 » 1 ,K .G.W*REL 

THROUGH ALICIA, FOR 1=1,1,1. G.M 

B(K,I )=X(K+1 , I) 

WHENEVER K.E.W, TRANSFER TO ALICIA 

A(K,I ]«UU + liI) 
ALICIA CONTINUF 

THROUGH MARTA, FOR K=l , 1 ,< .G. W*REL 

THROUGH MARTA, FOR I=1»1»I.G.M 

X(K,I)=R(K.I ) 

WHENEVER K.E.W, TRANSFER TO MART' 

U(K,I )=A(K,I > 
MARTA CONTINUE 

TRANSFER TO SONIA 
SILVIA PRINT COMMENT SEND OF EXECUTIONS 

PRINT COMMENT STO CONTINUE, GO TO THE TOP OF A NEW PAGES 

PRINT COMMENT SAND PRINT AN ASTERISKS 

READ DATA 

TRANSFER TO MAGDA 

END OF PROGRAM 
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DELFOR 

Purpose: to compute the plant transition matrices and the 
control transition matrices. 

Remarks: subroutine called by TIMDEL. 
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DELIA 
YOLIS 
MAGUE 



ROSANA 



MARIA 



SALOME 
ISABEL 
FANNY 



EXTERNAL FUNCTION (T) 

PROGRAM COMMON EM.DELF »M .R » W»A» B »D1 .D2 »U 

DIMENSION C(llOOO.H) ,A(400,G> »BU00»G> iEMUOOO»H» .XXU00»G) 

DIMENSION TERMUOO«G)»NTERMUOO.G) .UUUOO.G) .D1U00.G) 

DIMENSION DELFUOOO.H) .D2(400,G) .U(400,E) 

INTEGER I .J.K.L.M.N.Y.Q.R.W 

VECTOR VALUES H=3»1.0,0 

VECTOR VALUES G=2.1»0 

VECTOR VALUES E=2.1»0 

ENTRY TO DELFOR. 

G!2)=M 

H(2)=M 

H(3)=M 

LINDA=0. 

ROSA=-l. 

THROUGH YOLIS. FOR 1=1.1.1. G.M 

THROUGH YOLIS. FOR J=1.1»J.G.M 

WHENEVER J.G.R. TRANSFER TO DELIA 

Dl( I.J>=D1< I.J)*T 

D2(I »J>=D2( I»J)*T 

A( I.J)=A( I.J)*T 

TERM! I.J)=A< I »J) 

B( I.J)=B< I.J)*T 

N = 

MAXH=0. 

MAXV=0. 

THROUGH MARIA. FOR 1 = 1.1. 1. G.M 

SUMH=0. 

SUMV=0. 

THROUGH ROSANA. FOR J=1.1,J.G.M 

SUMH=SUMH+.ABS.TERM( I.J) 

SUMV=SUMV+.ABS.TERM< J.I ) 

CONTINUE 

WHENEVER SUMH.G.MAXH.MAXH=SUMH 

WHENEVER SUMV.G.MAXV.MAXV=SUMV 

CONTINUE 

NORM=MAXH 

WHENEVER MAXV.L.NORM ,NORM=MAXV 

WHENEVER LINDA. NE.O., TRANSFER TO CARMEN 

WHENEVER N.NE.O, TRANSFER TO CHELA 

SOLO=NORM 

K=2.*NORM 

WHENEVER K.L.2. 

W*l 

lN=K/2 

THROUGH SALOME. 



K=2 



FOR 1 = 1.1, I. G.M 



THROUGH SALOME. FOR J=1*1.J.G.M 

C(1.I.J)*0. 

WHENEVER J.E. I »C ( 1 . 1 » J > = 1. 

EM(W.I.J>*C(1.I»J> 

XX(I.J)=EM(W,I.J) 

UU(I.J)=0. 

TERM(I.J)=C(1»I.J) 

N=0 

WHENEVER N.GE.K. AND. LI NDA.E.O. , TRANSFER TO MAGUE 

WHENEVER N.GE.K. AND. LINDA. NE.O. .TRANSFER TO ELENA 

N*N+1 

THROUGH HILDA, FOR L=1»1»L-G.M 

THROUGH HILDA. FOR 1 = 1. 1,1. G.M 
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NTERMtL.I ) = 0. 

THROUGH LILIAf FOR J=1»1»J.G.M 

WHENEVER LINDA.F.0..C(N+1.J.I)=0. 

NTERMtL. I )=NTERM<L»I ) +A( L , J ) *TERM( J . I ) +B< L , J )*C (N+l »J» I ) 
t.ILIA CONTINUE 

EM1W.L.I >=EM( W.L.I )+NTERM(L»I )/(N+LINDA) 

WHENEVER R.E.O. TRANSFER TO HILDA 

XXlLiI >=XX(Li I)+NTERM(L»I)/( ( N+L INDA )* ( N+LINDA+1. ) > 
HILDA CONTINUE 

THROUGH JULIA. FOR 1 = 1.1. I. G.M 

THROUGH JULIA. FOR J=1»1.J.G.M 

CfN+l.I, J)=NTERM( I.J)/ (N+L INDA) 
JULIA TERMt I»J)=C(N+1»I.J> 

TRANSFER TO FANNY 
CHELA EPS=SOLO/(K+2.) 

RIJ=NORM*SOLO/( (K+1.)*(1.-EPS) ) 

THROUGH EVA. FOR 1=1*1. I. G.M 

THROUGH EVA. FOR J=1»1»J.G.M 

WW=.ABS. (EMfW.I ,J)*10..P.-07) 

WHENFVER RIJ.G.WW 

K=K+IN 

TRANSFER TO FANNY 

OTHERWISE 

TRANSFER TO EVA 

END OF CONDITIONAL 
EVA CONTINUF 
FLFNA LINDA=LINDA+1. 

ROSA=ROSA+l. 

WHENEVER R.E.O, TRANSFER TO PATY 

THROUGH MARTA. FOR L=1»1.L.G.M 

THROUGH MARTA. FOR 1=1.1. 1. G.R 

TERMIL.I >=0. 

THROUGH AURORA. FOR J=1.1»J.G.M 

TERMIL.I ) -TERM (L.I ) +XX < L . J ) *D1 1 J . I »+UU(Lt J>*D2 ( Jt H 
AURORA CONTINUE 
MARTA CONTINUE 

THROUGH IRMA. FOR 1=1.1.1. G.M 

THROUGH IRMA. FOR J=1.1.J.G.R 
IRMA DFLFI W.I »J)=TERM( I »J> 
PATY THROUGH SONIA. FOR 1 = 1.1. 1. G.M 

THROUGH SONIA. FOR J=1.1.J.G.M 

TERM! I,J>=EM<W,I »J> 
SONIA UIM I.JI=XX(I»J)' 

TRANSFER TO MAGUE 
CARMEN WHENEVER NORM. LE . 10 . .P. -07 .TRANSFER TO DIANA 

THROUGH YOCO. FOR L*ltltL*G.M 

THROUGH YOCO. FOR 1=1.1.1. G.M 

NTERM(L. I)=0. 

THROUGH JOSFFA. FOR J=1.1.J.G.M 
JOSE FA NTERM(L, I > = NTERM(L.I ) +B ( L . J )*C< 1 . J. I > 
YOCO CONTINUF 

W = W+1 

THROUGH ELISA. FOR 1=1.1.1. G.M 

THROUGH ELISA, FOR J=1.1.J.G.M 



C< 1 . I »J)=NTERM( I .J) /(ROSA+1. ) 
XXr>I»J)=C<l.I »J)/(ROSA + 2. ) 



EM(W,I.J)=C( l.I.J) 

ELISA TFRM(I,J)=C(1.I,J) 

TRANSFER TO ISABEL 
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PERTUR 

Purpose: to compute the forcing signal vector at the current 
time. The program keeps track of the past. 

Remarks: subroutine called by TIMDEL. 
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EXTERNAL FUNCTION (TA.LL) 

PROGRAM COMMON ' EM»DELFtMtRtWt.A»B»Ol*0.2 tU 

DIMENSION EM(4000,H> • DELFUOOO* H ) »A ( 400 »G) »B<40Q»G> 

D I MENS I ON 01 (AOO »G > »D2 ( 400 ,G > . U ( 400 » E ) 

INTEGER I*LL»R»W.M 

VECTOR VALUES G=2»l*0 

VECTOR VALUES E=2»l>0 

VECTOR VALUES H=3»ltOtO 

ENTRY TO PERTUR. 

G(2)=M 

H(2)=M 

E(2)=W 

H<3)=M 

U(LL»1)= 

U<LL»2)=" 

U(LL.R)= 

FUNCTION RETURN 
END OF FUNCTION 
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